{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 单纯形上的基函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 区间 $[x_0, x_1]$ 上 $p$ 次拉格朗日基函数的构造"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 重心坐标 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定区间单元上的一个重心坐标 $(\\lambda_0, \\lambda_1)$, 存在 $x \\in [x_0, x_1]$,  使得：\n",
    "\n",
    "$$\n",
    "\\lambda_0 := \\frac{x_1 - x}{x_1 - x_0}, \\quad\n",
    "\\lambda_1 := \\frac{x  - x_0}{x_1 - x_0}\n",
    "$$\n",
    "\n",
    "显然 \n",
    "$$ \n",
    "\\lambda_0 + \\lambda_1 = 1\n",
    "$$\n",
    "\n",
    "重心坐标关于  $x$ 的导数为：\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm d \\lambda_0}{\\mathrm d x} = -\\frac{1}{x_1 - x_0},\\quad \n",
    "\\frac{\\mathrm d \\lambda_1}{\\mathrm d x} = \\frac{1}{x_1 - x_0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "区间 $[x_0, x_1]$ 上的个 $p\\geq 1$ 次基函数共有 \n",
    "$$\n",
    "n_{dof} = p+1, \n",
    "$$ \n",
    "\n",
    "其计算公式如下:\n",
    "\n",
    "$$\n",
    "\\phi_{m,n} = \\frac{p^p}{m!n!}\\prod_{l_0 = 0}^{m - 1}\n",
    "(\\lambda_0 - \\frac{l_0}{p}) \\prod_{l_1 = 0}^{n-1}(\\lambda_1 -\n",
    "\\frac{l_1}{p}).\n",
    "$$\n",
    "\n",
    "其中 $ m\\geq 0$, $n\\geq 0$, 且 $m+n=p$, 这里规定:\n",
    "\n",
    "$$\n",
    " \\prod_{l_i=0}^{-1}(\\lambda_i - \\frac{l_i}{p}) := 1,\\quad i=0, 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $p$ 次基函数的面向数组的计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造向量：\n",
    "$$\n",
    "P = ( \\frac{1}{0!},  \\frac{1}{1!}, \\frac{1}{2!}, \\cdots, \\frac{1}{p!})\n",
    "$$\n",
    "\n",
    "构造矩阵：\n",
    "$$\n",
    "A :=                                                                            \n",
    "\\begin{pmatrix}  \n",
    "1  &  1  \\\\\n",
    "\\lambda_0 & \\lambda_1\\\\                                             \n",
    "\\lambda_0 - \\frac{1}{p} & \\lambda_1 - \\frac{1}{p}\\\\   \n",
    "\\vdots & \\vdots \\\\                                                     \n",
    "\\lambda_0 - \\frac{p - 1}{p} & \\lambda_1 - \\frac{p - 1}{p}\n",
    "\\end{pmatrix}                                                                   \n",
    "$$ \n",
    "\n",
    "对 $A$ 的每一列做累乘运算, 并左乘由 $P$ 形成的对角矩阵, 得矩阵:\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1\\\\\n",
    "\\lambda_0 & \\lambda_1\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "易知, 只需从 $B$ 的每一列中各选择一项相乘(要求二项次数之和为 $p$), 再乘以 $p^p$ 即可得到相应的基函数, 其中取法共有 \n",
    "\n",
    "$$\n",
    "n_{dof} = {p+1}\n",
    "$$\n",
    "\n",
    "构造指标矩阵：\n",
    "\n",
    "$$\n",
    "I = \\begin{pmatrix}\n",
    "p  & 0 \\\\ p-1 & 1 \\\\ \\vdots & \\vdots \\\\ 0 & p \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "则 $p+1$ 个 $p$ 次基函数可写成如下形式 \n",
    "\n",
    "$$\n",
    "\\phi_i = p^pB_{I_{i,0}, 0}B_{I_{i, 1},1}, \\quad i = 0, 1, \\cdots, n_{dof}\n",
    "$$\n",
    "\n",
    "对应 `Python` 的代码如下:\n",
    "\n",
    "```Python\n",
    "phi = p**p*np.prod(B[I, [0, 1]], axis=1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造矩阵:\n",
    "$$\n",
    "F^i = \n",
    "\\begin{pmatrix}\n",
    "1 & \\lambda_i & \\cdots & \\lambda_i \\\\\n",
    "\\lambda_i - \\frac{1}{p} & 1 & \\cdots & \\lambda_i - \\frac{1}{p} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\lambda_i - \\frac{p-1}{p} & \\lambda_i - \\frac{p-1}{p} & \\cdots & 1\n",
    "\\end{pmatrix}_{p\\times p}\n",
    ", \\quad i = 0, 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记\n",
    "$$\n",
    "f_{i,j} = \\sum_{m=0}^j\\prod_{k=0}^j F^i_{k, m},\\quad i = 0, 1, 0 \\leq j \\leq p-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla B = &  \n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 \\\\\n",
    "\\nabla \\lambda_1 & \\nabla \\lambda_2 \\\\\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p}) & \n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p})  \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 \\\\\n",
    "\\nabla \\lambda_1 & \\nabla \\lambda_2\\\\\n",
    "f_{0,0} \\nabla \\lambda_1 & \n",
    "f_{1,0} \\nabla \\lambda_2  \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "f_{0, p-1} \\nabla \\lambda_1 &\n",
    "f_{1, p-1} \\nabla \\lambda_2 \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0\\\\\n",
    "1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} \\\\\n",
    "\\vdots & \\vdots  \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1}   \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\nabla \\lambda_0 & 0 \\\\\n",
    "0 & \\nabla \\lambda_1 \n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 三角形单元"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 重心坐标"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定三角形单元 $\\tau$, 其三个顶点 $\\mathbf x_0 :=(x_0,y_0)$, $\\mathbf x_1 :=(x_1,y_1)$ 和 $\\mathbf x_2 :=(x_2,y_2)$ 逆时针排列, 且不在同一条直线上, 那么向量 $\\vec{\\mathbf x_0\\mathbf x_1}$ 和 $\\vec{\\mathbf x_0\\mathbf x_2}$ 是线性无关的. 这等价于矩阵\n",
    "\n",
    "$$\n",
    "A = \n",
    "\\begin{pmatrix}\n",
    "x_0 & x_1 & x_2 \\\\\n",
    "y_0 & y_1 & y_2 \\\\\n",
    "1   & 1   & 1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "非奇异. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "任给一点 $\\mathbf{x}:=(x,y)\\in\\tau$, 求解下面的线性方程组\n",
    "\n",
    "$$\n",
    "A \n",
    "\\begin{pmatrix}\n",
    "\\lambda_0 \\\\\n",
    "\\lambda_1\\\\\n",
    "\\lambda_2  \n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "x \\\\\n",
    "y\\\\\n",
    "1  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "可得唯一的一组解$\\lambda_0,\\lambda_1,\\lambda_2$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此对任意二维点 $\\mathbf{x}\\in\\tau$, 有\n",
    "\n",
    "$$\n",
    "\\mathbf{x}=\\lambda_0(\\mathbf{x})\\mathbf{x}_0 + \\lambda_1(\\mathbf{x})\\mathbf{x}_1 + \\lambda_2(\\mathbf{x})\\mathbf{x}_2 \n",
    "\\text{ 且 } \\lambda_0(\\mathbf{x}) + \\lambda_1(\\mathbf{x}) + \\lambda_2(\\mathbf{x}) = 1. \n",
    "$$\n",
    "\n",
    "$\\lambda_0,\\lambda_1,\\lambda_2$ 称为点 $\\mathbf{x}$ 关于点 $\\mathbf{x}_0,\\mathbf{x}_1$ 和$\\mathbf{x}_2$ 的**重心坐标**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, $\\lambda_0, \\lambda_1, \\lambda_2$ 都是关于 $\\mathbf x$ 的线性函数, 且有\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\lambda_0(\\mathbf x_0) = 1,\\quad & \\lambda_0(\\mathbf x_1) = 0,\\quad& \\lambda_0(\\mathbf x_2) = 0\\\\\n",
    "\\lambda_1(\\mathbf x_0) = 0,\\quad & \\lambda_1(\\mathbf x_1) = 1,\\quad& \\lambda_1(\\mathbf x_2) = 0\\\\\n",
    "\\lambda_2(\\mathbf x_0) = 0,\\quad & \\lambda_2(\\mathbf x_1) = 0,\\quad & \\lambda_2(\\mathbf x_2) = 1\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\lambda_0, \\lambda_1, \\lambda_2$ 关于 $\\mathbf x$ 的梯度分别为:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla\\lambda_0 = \\frac{1}{2|\\tau|}(\\mathbf x_2 - \\mathbf x_1)W\\\\\n",
    "\\nabla\\lambda_1 = \\frac{1}{2|\\tau|}(\\mathbf x_0 - \\mathbf x_2)W\\\\\n",
    "\\nabla\\lambda_2 = \\frac{1}{2|\\tau|}(\\mathbf x_1 - \\mathbf x_0)W\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "其中 \n",
    "\n",
    "$$\n",
    "W = \\begin{pmatrix}\n",
    "0 & 1\\\\ -1 & 0 \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  三角形单元上的 $p$ 次基函数公式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定三角形单元上的一个重心坐标 $(\\lambda_0, \\lambda_1, \\lambda_2)$, 所有 $p\\geq 1$ 次基函数的计算公式如下:\n",
    "\n",
    "$$\n",
    "\\phi_{m,n,k} = \\frac{p^p}{m!n!k!}\\prod_{l_0 = 0}^{m - 1}\n",
    "(\\lambda_0 - \\frac{l_0}{p}) \\prod_{l_1 = 0}^{n-1}(\\lambda_1 -\n",
    "\\frac{l_1}{p}) \\prod_{l_2=0}^{k-1}(\\lambda_2 - \\frac{l_2}{p}).\n",
    "$$\n",
    "\n",
    "其中 $ m\\geq 0$, $n\\geq 0$, $ k \\geq 0$, 且 $m+n+k=p$, 这里规定:\n",
    "\n",
    "$$\n",
    " \\prod_{l_i=0}^{-1}(\\lambda_i - \\frac{l_i}{p}) := 1,\\quad i=0, 1, 2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 三角形单元上的 $p$ 次基函数面向矩阵计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造向量：\n",
    "$$\n",
    "P = ( \\frac{1}{0!},  \\frac{1}{1!}, \\frac{1}{2!}, \\cdots, \\frac{1}{p!})\n",
    "$$\n",
    "\n",
    "构造矩阵：\n",
    "$$\n",
    "A :=                                                                            \n",
    "\\begin{pmatrix}  \n",
    "1  &  1  & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\                                             \n",
    "\\lambda_0 - \\frac{1}{p} & \\lambda_1 - \\frac{1}{p} & \\lambda_2 - \\frac{1}{p}\\\\   \n",
    "\\vdots & \\vdots & \\vdots \\\\                                                     \n",
    "\\lambda_0 - \\frac{p - 1}{p} & \\lambda_1 - \\frac{p - 1}{p} & \\lambda_2 - \\frac{p - 1}{p}\n",
    "\\end{pmatrix}                                                                   \n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 $A$ 的每一列做累乘运算, 并左乘由 $P$ 形成的对角矩阵, 得矩阵:\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, 只需从 $B$ 的每一列中各选择一项相乘(要求三项次数之和为 $p$), 再乘以 $p^p$ 即可得到相应的基函数, 其中取法共有 \n",
    "\n",
    "$$\n",
    "n_{dof} = {(p+1)(p+2)\\over 2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "即需要构造 $n_{dof}\\times 3$ 的指标矩阵 $I$, \n",
    "\n",
    "$$\n",
    "I = \\begin{pmatrix}\n",
    "p & 0 & 0 \\\\\n",
    "p-1 & 1 & 0 \\\\\n",
    "p-1 & 0 & 1 \\\\\n",
    "p-2 & 2 & 0 \\\\\n",
    "p-2 & 1 & 1 \\\\\n",
    "p-2 & 0 & 2 \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "0   & p & 0 \\\\\n",
    "0   & p-1 & 1\\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "0   & 0  & p\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "则基函数可写成如下形式 \n",
    "\n",
    "$$\n",
    "\\phi_i = p^pB_{I_{i,0}, 0}B_{I_{i, 1},1}B_{I_{i, 2}, 2}, \\quad i = 0, 1, \\cdots, n_{dof}\n",
    "$$\n",
    "\n",
    "对应 `Python` 的代码如下:\n",
    "\n",
    "```Python\n",
    "phi = p**p*np.prod(B[I, [0, 1, 2]], axis=1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 三角形单元上的 $p$ 次基函数梯度面向数组的计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**分析:**\n",
    "\n",
    "根据求导法则, 关键是计算出矩阵\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "中每个元素的梯度. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造矩阵:\n",
    "$$\n",
    "F^i = \n",
    "\\begin{pmatrix}\n",
    "1 & \\lambda_i & \\cdots & \\lambda_i \\\\\n",
    "\\lambda_i - \\frac{1}{p} & 1 & \\cdots & \\lambda_i - \\frac{1}{p} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\lambda_i - \\frac{p-1}{p} & \\lambda_i - \\frac{p-1}{p} & \\cdots & 1\n",
    "\\end{pmatrix}_{p\\times p}\n",
    ", \\quad i = 0, 1, 2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记\n",
    "$$\n",
    "f_{i,j} = \\sum_{m=0}^j\\prod_{k=0}^j F^i_{k, m},\\quad i = 0, 1, 2, 0 \\leq j \\leq p-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla B = &  \n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "\\nabla \\lambda_0 & \\nabla \\lambda_1 & \\nabla \\lambda_2\\\\\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0 \\\\\n",
    "\\nabla \\lambda_0 & \\nabla \\lambda_1 & \\nabla \\lambda_2\\\\\n",
    "f_{0,0} \\nabla \\lambda_0 & \n",
    "f_{1,0} \\nabla \\lambda_1 &\n",
    "f_{2,0} \\nabla \\lambda_2 \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} \\nabla \\lambda_0 &\n",
    "f_{1, p-1} \\nabla \\lambda_1 &\n",
    "f_{2, p-1} \\nabla \\lambda_2 \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1}  \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\nabla \\lambda_0 & 0 & 0\\\\\n",
    "0 & \\nabla \\lambda_1 & 0 \\\\\n",
    "0 & 0 & \\nabla \\lambda_2\n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记 \n",
    "$$\n",
    "F =\n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1}  \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用前面的矩阵 $B$, $I$, $Dlambda$, 如下代码即可计算所有单元上的 $p$ 次元基函数的导数:\n",
    "\n",
    "```Python\n",
    "pp = p**p\n",
    "ldof = (p+1)*(p+2)/2\n",
    "Q = B[I,[0,1,2]]\n",
    "M = F[I,[0,1,2]]\n",
    "R = np.zeros((ldof,3), dtype=np.float)\n",
    "R[:,0] = M[:,0]*Q[:,1]*Q[:,2]\n",
    "R[:,1] = Q[:,0]*M[:,1]*Q[:,2]\n",
    "R[:,2] = Q[:,0]*Q[:,1]*M[:,2]\n",
    "gradphi = np.einsum('ij, kj...->ki...', pp*R, Dlambda) # 爱因斯坦求和\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 四面体单元"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 重心坐标"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定四面体单元 $\\tau$, 其三个顶点 $\\mathbf x_0 :=(x_0, y_0，z_0)$, $\\mathbf x_1 :=(x_1,y_1,z_1)$, $\\mathbf x_2 :=(x_2,y_2,z_2)$ 和 $\\mathbf x_3 :=(x_3, y_3, z_3)$ 排列满足右手法则, 体积不为零. 这等价于矩阵\n",
    "\n",
    "$$\n",
    "A = \n",
    "\\begin{pmatrix}\n",
    "x_0 & x_1 & x_2 & x_3\\\\\n",
    "y_0 & y_1 & y_2 & y_3\\\\\n",
    "z_0 & x_1 & x_2 & x_3\\\\\n",
    "1   & 1   & 1 & 1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "非奇异. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "任给一点 $\\mathbf{x}:=(x,y)\\in\\tau$, 求解下面的线性方程组\n",
    "\n",
    "$$\n",
    "A \n",
    "\\begin{pmatrix}\n",
    "\\lambda_0 \\\\\n",
    "\\lambda_1\\\\\n",
    "\\lambda_2 \\\\\n",
    "\\lambda_3\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "x \\\\\n",
    "y\\\\\n",
    "z \\\\\n",
    "1  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "可得唯一的一组解$\\lambda_0,\\lambda_1,\\lambda_2, \\lambda_3$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此对任意三维点 $\\mathbf{x}\\in\\tau$, 有\n",
    "\n",
    "$$\n",
    "\\mathbf{x}=\\lambda_0(\\mathbf{x})\\mathbf{x}_0 + \\lambda_1(\\mathbf{x})\\mathbf{x}_1 + \\lambda_2(\\mathbf{x})\\mathbf{x}_2 + \\lambda_3(\\mathbf{x})\\mathbf{x}_3 \n",
    "\\text{ 且 } \\lambda_0(\\mathbf{x}) + \\lambda_1(\\mathbf{x}) + \\lambda_2(\\mathbf{x}) + \\lambda_3(\\mathbf{x}) = 1. \n",
    "$$\n",
    "\n",
    "$\\lambda_0,\\lambda_1,\\lambda_2, \\lambda_3$ 称为点 $\\mathbf{x}$ 关于点 $\\mathbf{x}_0,\\mathbf{x}_1$, $\\mathbf{x}_2$ 和 $\\mathbf x_3$ 的**重心坐标**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, $\\lambda_0, \\lambda_1, \\lambda_2， \\lambda_3$ 都是关于 $\\mathbf x$ 的线性函数, 且有\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\lambda_0(\\mathbf x_0) = 1,\\quad & \\lambda_0(\\mathbf x_1) = 0,\\quad& \\lambda_0(\\mathbf x_2) = 0,\\quad & \\lambda_0(\\mathbf x_3) = 0\\\\\n",
    "\\lambda_1(\\mathbf x_0) = 0,\\quad & \\lambda_1(\\mathbf x_1) = 1,\\quad& \\lambda_1(\\mathbf x_2) = 0,\\quad & \\lambda_0(\\mathbf x_3) = 0\\\\\n",
    "\\lambda_2(\\mathbf x_0) = 0,\\quad & \\lambda_2(\\mathbf x_1) = 0,\\quad & \\lambda_2(\\mathbf x_2) = 1,\\quad & \\lambda_0(\\mathbf x_3) = 0\\\\\n",
    "\\lambda_3(\\mathbf x_0) = 0,\\quad & \\lambda_3(\\mathbf x_1) = 0,\\quad & \\lambda_3(\\mathbf x_2) = 0,\\quad & \\lambda_0(\\mathbf x_3) = 1\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "三维情形的基函数表示如下:\n",
    "$$\n",
    "\\phi_{m,n,k,v} = \\frac{p^p}{m!n!k!v!}\\prod_{l_0 = 0}^{m - 1}\n",
    "(\\lambda_0 - \\frac{l_0}{p}) \\prod_{l_1 = 0}^{n-1}(\\lambda_1 -\n",
    "\\frac{l_1}{p}) \\prod_{l_2=0}^{k-1}(\\lambda_2 - \\frac{l_2}{p})\\prod_{l_3=0}^{v-1}(\\lambda_3 - \\frac{l_3}{p}).\n",
    "$$\n",
    "\n",
    "其中 $ m\\geq 0$, $n\\geq 0$, $ k \\geq 0$, $v\\geq 0$,  且 $m+n+k+v=p$, 这里规定:\n",
    "\n",
    "$$\n",
    " \\prod_{l_i=0}^{-1}(\\lambda_i - \\frac{l_i}{p}) := 1,\\quad i=0, 1, 2, 3\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造向量：\n",
    "$$\n",
    "P = ( \\frac{1}{0!},  \\frac{1}{1!}, \\frac{1}{2!}, \\cdots, \\frac{1}{p!})\n",
    "$$\n",
    "\n",
    "构造矩阵：\n",
    "$$\n",
    "A :=                                                                            \n",
    "\\begin{pmatrix}  \n",
    "1  &  1  & 1 & 1\\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2 & \\lambda_3\\\\                                             \n",
    "\\lambda_0 - \\frac{1}{p} & \\lambda_1 - \\frac{1}{p} & \\lambda_2 - \\frac{1}{p} & \\lambda_3 - \\frac{1}{p}\\\\   \n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\                                                     \n",
    "\\lambda_0 - \\frac{p - 1}{p} & \\lambda_1 - \\frac{p - 1}{p} & \\lambda_2 - \\frac{p - 1}{p}& \\lambda_3 - \\frac{p - 1}{p}\n",
    "\\end{pmatrix}                                                                   \n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 $A$ 的每一列做累乘运算, 并左乘由 $P$ 形成的对角矩阵, 得矩阵:\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 & 1\\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2 & \\lambda_3\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_3 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_3 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, 只需从 $B$ 的每一列中各选择一项相乘(要求四项次数之和为 $p$), 再乘以 $p^p$ 即可得到相应的基函数, 其中取法共有 \n",
    "\n",
    "$$\n",
    "n_{dof} = {(p+1)(p+2)(p+3)\\over 6}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "即需要构造 $n_{dof}\\times 3$ 的指标矩阵 $I$, \n",
    "\n",
    "$$\n",
    "I = \\begin{pmatrix}\n",
    "p & 0 & 0 & 0\\\\\n",
    "p-1 & 1 & 0 & 0\\\\\n",
    "p-1 & 0 & 1 & 0\\\\\n",
    "p-1 & 0 & 0 & 1\\\\\n",
    "p-2 & 2 & 0 & 0\\\\\n",
    "p-2 & 1 & 1 & 0\\\\\n",
    "p-2 & 1 & 0 & 1\\\\\n",
    "p-2 & 0 & 2 & 0 \\\\\n",
    "p-2 & 0 & 1 & 1 \\\\\n",
    "p-2 & 0 & 0 & 2\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "0   & p & 0 & 0 \\\\\n",
    "0   & p-1 & 1 & 0\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "0   & 0  & 0 & p \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "则基函数可写成如下形式 \n",
    "\n",
    "$$\n",
    "\\phi_i = p^pB_{I_{i,0}, 0}B_{I_{i, 1},1}B_{I_{i, 2}, 2}B_{I_{i, 3}, 3}, \\quad i = 0, 1, \\cdots, n_{dof}\n",
    "$$\n",
    "\n",
    "对应 `Python` 的代码如下:\n",
    "\n",
    "```Python\n",
    "phi = p**p*np.prod(B[I, [0, 1, 2, 3]], axis=1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 四面体单元上的 $p$ 次基函数梯度面向数组的计算"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**分析:**\n",
    "\n",
    "根据求导法则, 关键是计算出矩阵\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 & 1\\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2& \\lambda_3\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p})& \\prod_{l=0}^{1}(\\lambda_3 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p})& \\prod_{l=0}^{p-1}(\\lambda_3 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "中每个元素的梯度. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造矩阵:\n",
    "$$\n",
    "F^i = \n",
    "\\begin{pmatrix}\n",
    "1 & \\lambda_i & \\cdots & \\lambda_i \\\\\n",
    "\\lambda_i - \\frac{1}{p} & 1 & \\cdots & \\lambda_i - \\frac{1}{p} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\lambda_i - \\frac{p-1}{p} & \\lambda_i - \\frac{p-1}{p} & \\cdots & 1\n",
    "\\end{pmatrix}_{p\\times p}\n",
    ", \\quad i = 0, 1, 2, 3\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记\n",
    "$$\n",
    "f_{i,j} = \\sum_{m=0}^j\\prod_{k=0}^j F^i_{k, m},\\quad i = 0, 1, 2,3\\quad 0 \\leq j \\leq p-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla B = &  \n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 &0 & 0 & 0\\\\\n",
    "\\nabla \\lambda_0 &\\nabla \\lambda_1 & \\nabla \\lambda_2 & \\nabla \\lambda_3\\\\\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p}) & \n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_3 - \\frac{l}{p}) \\\\\n",
    "\\vdots &\\vdots & \\vdots & \\vdots \\\\\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_3 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 &0 & 0 & 0 \\\\\n",
    "\\nabla \\lambda_0 &\\nabla \\lambda_1 & \\nabla \\lambda_2 & \\nabla \\lambda_3\\\\\n",
    "f_{0,0} \\nabla \\lambda_0 & \n",
    "f_{1,0} \\nabla \\lambda_1 &\n",
    "f_{2,0} \\nabla \\lambda_2 &\n",
    "f_{3,0} \\nabla \\lambda_3 \\\\\n",
    "\\vdots &\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} \\nabla \\lambda_0 &\n",
    "f_{1, p-1} \\nabla \\lambda_1 &\n",
    "f_{2, p-1} \\nabla \\lambda_2 &\n",
    "f_{3, p-1} \\nabla \\lambda_3 \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 &0 & 0 & 0\\\\\n",
    "1 &1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} &\n",
    "f_{3,0} \\\\\n",
    "\\vdots &\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1} &\n",
    "f_{3, p-1}\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\nabla \\lambda_0 & 0 & 0& 0\\\\\n",
    "0 & \\nabla \\lambda_1 & 0 & 0\\\\\n",
    "0 & 0 & \\nabla \\lambda_2 & 0\\\\\n",
    "0 & 0 & 0 & \\nabla \\lambda_3 \n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记 \n",
    "$$\n",
    "F =\n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0& 0\\\\\n",
    "1 & 1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} &\n",
    "f_{3,0} \\\\\n",
    "\\vdots &\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1} &\n",
    "f_{3, p-1}  \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用前面的矩阵 $B$, $I$, $Dlambda$, 如下代码即可计算所有单元上的 $p$ 次元基函数的导数:\n",
    "\n",
    "```Python\n",
    "pp = p**p\n",
    "ldof = p(p+1)*(p+2)/6\n",
    "Q = B[I,[0,1,2,3]]\n",
    "M = F[I,[0,1,2,3]]\n",
    "R = np.zeros((ldof,4), dtype=np.float)\n",
    "R[:,0] = M[:,0]*Q[:,1]*Q[:,2]*Q[:,3]\n",
    "R[:,1] = Q[:,0]*M[:,1]*Q[:,2]*Q[:,3]\n",
    "R[:,2] = Q[:,0]*Q[:,1]*M[:,2]*Q[:,3]\n",
    "R[:,3] = Q[:,0]*Q[:,1]*Q[:,2]*M[:,3]\n",
    "gradphi = np.einsum('ij, kj...->ki...', pp*R, Dlambda) # 爱因斯坦求和\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 四边形单元"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定一个四边形单元 $\\tau := (\\mathbf x_0, \\mathbf x_1, \\mathbf x_2, \\mathbf x_3)$, 可以建立它与参考单元之间的映射关系\n",
    "\n",
    "$$\n",
    "A :=\n",
    "\\begin{pmatrix}\n",
    "1 - \\xi \\\\\n",
    "\\xi\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "1 - \\eta & \\eta\n",
    "\\end{pmatrix}\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "(1-\\xi)(1-\\eta) & (1-\\xi)\\eta\\\\\n",
    "\\xi(1-\\eta) & \\xi\\eta\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\mathbf x(\\xi, \\eta) =\n",
    "\\begin{pmatrix}\n",
    "(1-\\xi)(1-\\eta) & (1-\\xi)\\eta & \\xi(1-\\eta) & \n",
    "\\xi\\eta\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "x_0 \\\\ x_1 \\\\ x_2 \\\\ x_3\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\frac{\\mathrm d x}{\\mathrm d \\xi} & \\frac{\\mathrm d x}{\\mathrm d\\eta}\\\\\n",
    "\\frac{\\mathrm d y}{\\mathrm d \\xi} & \\frac{\\mathrm d y}{\\mathrm d\\eta}\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "(\\eta - 1)x_0 -\\eta x_1 + (1 - \\eta) x_2 + \\eta x_3 & (\\xi - 1)x_0 + (1 - \\xi) x_1  - \\xi x_2 + \\xi x_3\\\\\n",
    "(\\eta - 1)y_0 -\\eta y_1 + (1 - \\eta) y_2 + \\eta y_3 & (\\xi - 1)y_0 + (1 - \\xi) y_1  - \\xi y_2 + \\xi y_3\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\frac{\\mathrm d x}{\\mathrm d \\xi} & \\frac{\\mathrm d x}{\\mathrm d\\eta}\\\\\n",
    "\\frac{\\mathrm d y}{\\mathrm d \\xi} & \\frac{\\mathrm d y}{\\mathrm d\\eta}\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "x_2 - x_0 & x_1 - x_0\\\\\n",
    "y_2 - y_0 & y_1 - y_0\n",
    "\\end{pmatrix}\n",
    "+\n",
    "\\begin{pmatrix}\n",
    "x_0 - x_1 - x_2 + x_3\\\\\n",
    "y_0 - y_1 - y_2 + y_3\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\eta & \\xi\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记 $\\lambda_0 = 1-x$, $\\lambda_1 = x$, 利用上面的方法可以构造 $[0, 1]$ 区间上的一元 $p$ 次元基函数\n",
    "\n",
    "\n",
    "对应 `Python` 的代码如下:\n",
    "\n",
    "```Python\n",
    "phi = p**p*np.prod(B[I, [0, 1]], axis=1)\n",
    "```\n",
    "\n",
    "每个基函数关于 $x$ 的导数的计算代码如下：\n",
    "\n",
    "\n",
    "```Python\n",
    "pp = p**p\n",
    "ldof = p+1\n",
    "Q = B[I,[0,1]]\n",
    "M = F[I,[0,1]]\n",
    "gradphi = -M[:,0]*Q[:,1] + Q[:,0]*M[:,1]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$$\n",
    "\\Phi(x) = (\\phi_0(x), \\phi_1(x), \\cdots, \\phi_{p}(x))\n",
    "$$\n",
    "\n",
    "参考单元上的双 $p$ 次元基函数为 \n",
    "\n",
    "$$\n",
    "\\Phi^T(\\xi)\\Phi(\\eta) =\n",
    "\\begin{pmatrix}\n",
    "\\phi_0(\\xi)\\phi_0(\\eta) & \\phi_0(\\xi)\\phi_1(\\eta) & \\cdots & \\phi_0(\\xi)\\phi_p(\\eta)\\\\\n",
    "\\phi_1(\\xi)\\phi_0(\\eta) & \\phi_1(\\xi)\\phi_1(\\eta) & \\cdots & \\phi_1(\\xi)\\phi_p(\\eta)\\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\phi_p(\\xi)\\phi_0(\\eta) & \\phi_p(\\xi)\\phi_1(\\eta) & \\cdots & \\phi_p(\\xi)\\phi_p(\\eta)\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "梯度的计算公式为\n",
    "\n",
    "$$\n",
    "\\left(\\frac{\\mathrm d \\Phi^T(\\xi)}{\\mathrm d \\xi}\\Phi(\\eta), \\Phi^T(\\xi)\\frac{\\mathrm d \\Phi^T(\\eta)}{\\mathrm d\\eta}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 六面体单元"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定一个六面体单元 $\\tau := (x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7)$, 建立它与参考单元的映射关系："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUkAAAFKCAYAAACQHq0GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzsvXl4VPXZ//8+s2WdSWbJHrJAIAmBhLAFEBDECghE0afu\n8vh4tda2XqL10Yr2K1W8rHuLVltbFR4rQl0o1lZUloBAQGKEbGSbLJNJMklIJskkmSSTmTO/P/jN\naZYh62zn5H5dF39wZuacz5mZvOf+fO77874Zu90OgiAIwjkibw+AIAjClyGRJAiCGAUSSYIgiFEg\nkSQIghgFEkmCIIhRIJEkCIIYBRJJgiCIUSCRJAiCGAUSSYIgiFGQTOC5tDWHIAihwYz1BIokCYIg\nRoFEkiAIYhRIJAmCIEaBRJIgCGIUSCQJgiBGgUSSIAhiFEgkCYIgRoFEkiAIYhRIJAmCIEaBRJKY\nMiaTCU899RSSk5MRFxeHrKwsfP7559zjn376KdLT0xEfH4/MzEx8+OGHXhwtQUwMEkliStTV1WHl\nypXQaDQoKipCXV0dNmzYgFtuuQX79u3Du+++iz//+c/IyclBRUUF5s6di23btuHo0aPeHjpBjAtm\nAt0Sae82MYSBgQEsW7YMzz77LDZv3swdN5vNkMvl0Gg0CA4OxoULF6BQKPDYY4/h97//PRiGwauv\nvopHH33Ui6MnCAC0d5twJ2+88QaysrKGCCQAiMViAEBrayvuv/9+KBQKAIDBYADDMJDL5bjppptG\nnO+LL77AunXrsHr1asyfPx8PP/ww6uvr3X8jBDEKFEkSk2bWrFk4ceIEZsyYMeR4fn4+lixZAoZh\nUFRUhLlz5wIAOjs78fXXX2P58uUjXvOXv/wF+/fvx8cff4ywsDB0dnZi2bJlaGpqwsmTJ5Genu6x\n+yKmFWNGkrDb7eP9RxAcAwMD9n379jl97I033rAzDGOPiIgY17lMJpN9xowZdqPROOT4m2++aWcY\nxr527dopj5cgrsKY2kfTbWJSSCQS3HXXXU4fy8nJAQCsXr16XOc6c+YM6uvrR5xv9uzZAIDvvvtu\nCiMliKlBIkm4FLvdjpMnT4JhGKxZs2ZcrzGbzQCAb775BkajkTtusVgAgFvTJAhvMBFncoIYk4KC\nArS3t09IJDdu3IgtW7Zg1qxZUKlU3PGSkhIAwLp169wxVIIYFySShEtxTLU1Gg2XsBmMVquF0WjE\n0qVLuWMBAQFDis8dHDx4EAEBAfjNb37jvgETxBjQdJuYNEeOHEFeXt6QYw6RvOaaa5y+5umnn0ZZ\nWdmY53777bdRUVGBzz77DCkpKVMfLEFMEookiUmxfft2vPnmmxCLxSgtLUVSUhK6urpw7NgxMAyD\ntLS0Ea/RarXIzc3F3r17nZ7TZrNh48aNMBqN0Gq12L17NzZs2ODmOyGI0aFIkpgUe/fuBcMw8PPz\n4xIrr7/+OjQaDQDAarUOeX5/fz/uvvtuvPLKKwgICHB6TrFYjG+++Qbff/89ysrK8Nprr+Gaa65B\nY2Oje2+GIEaBRJKYFCkpKVAoFPjHP/6B8PBw7Nu3D4cOHUJeXh5SUlLw0UcfoaGhAQBQVlaGDRs2\nIDs7G3fccce4zh8ZGYkXX3wRZ8+excaNG7lMN0F4GtpxQ0yKhoYGPP7448jPz8fAwABWrVqF1157\nDRqNBkajETt27MDhw4chkUgQGxuLJ598EjfeeOOErmGxWBAUFASWZbFnzx5s27bNTXdDTGPG3HFD\nIkmMC4vFgsOHD6O9vR3z5s1DZGQkIiMjIZFMfVn7D3/4Az755BO8/PLLIxI+UVFRaGlpwS9/+Uu8\n8cYbU74WQQxjTJGkxA0xJi0tLTh+/Djy8vJgt9vR0NAAqVSK4OBgxMbGIjo6ekqi+Zvf/Aa9vb3c\nGuRgent7AQAhISEuuReCmCgkksSoFBcX4+TJk+jq6kJERAQAYMGCBejp6YHRaERlZSWKioogk8k4\n0YyJiUFkZCTCw8PHJZoxMTFobGzEbbfdNuS4wWCAyWQCwzCU5Sa8Bokk4ZS+vj7k5OSgsLAQcrkc\nCxYsQEVFBRfZBQUFISgoiHPz6erqQnt7O8rLy1FUVASpVAq5XI64uDgu0gwPD4dINDJXuH37dnz7\n7bf48Y9/POT4vn37AAC33377VesuCcLd0JokMYLGxkYcP34ctbW1mDVrFtRqNQCgvLwcZrMZmZmZ\nY57DIZodHR3o7++HTCaDQqFAXFwcoqKiRojmK6+8gn/84x+45557kJSUhNLSUjzzzDO499578frr\nr0Mmk7n1nolpCyVuiPHDsix++OEHnDlzBv39/UhJSYG/vz/3eHl5OXp6erBw4cIJn7e7uxtGoxEm\nkwn9/f3w9/eHXC7HjBkzuOm5WCzG0aNHodPpoFarccMNN4zwnSQIF0MiSYyP7u5u5OTkoKSkBEql\nEomJiSOeU1FRga6uLixatGhK12JZFiaTCe3t7ejs7MTAwAD8/f0REhKC+Ph4LtJUqVROp+cE4UJI\nJImx0el0OHr0KJqamkY48QymsrISnZ2dWLx4sUuvP1g0TSYTLBYLJ5qJiYmIjIxEVFTUVcdFEFOA\nSoCIq2O1WpGXl4ezZ8/CbrdjwYIFkEqlHh+HSCRCaGgoQkNDAVwRzY6ODnR0dCA/Px8WiwWBgYFQ\nKpWIi4sj0SQ8ConkNKWjowPHjh1DaWkpoqKifGrtTyQSQaVScSJos9nQ2dkJo9GI/Px8WK1W+Pv7\nQ6VSISEhAREREYiKiuJEliBcCYnkNKSyshI5OTlobW1Famoq5HL5uF7HMAwmsDzjMsRi8QjR7Ojo\nQHt7O86dOweWZREYGAi1Wo34+HhONMnRnHAFJJLTCIvFgtzcXOTl5UEmkyEzM5Nr/8onxGIx1Go1\nV5o0WDRzc3Nhs9kQFBQEjUaDuLg4Ek1iSpBIThNaW1tx9OhRVFVVISYmBtHR0RM+B8OM3X3TGwwX\nzYGBAU406+vrwbIsJ5qDI83g4GAvj5zgAySS04Di4mJ8++236OzsRFpaGgIDAyd1Hm9NtyeKVCpF\nWFgYwsLCAFwRzfb2dhiNRtTV1QEAAgMDERERgRkzZnCJoMm+L4SwIZEUMH19fTh16hR++OEHBAUF\nITMzc0p1h3wRyeFIpVKEh4cjPDwcwJVlh46ODly+fBm1tbWw2+0IDAxEZGQklz2PiIgg0SQAkEgK\nlsFbCxMTE7moaqrwUSSHI5PJRoim0WhEU1MTqqurAVzZmx4RETGk5Gjw7iNi+kAiKTBYlsXFixdx\n+vRp9PX1IT093WV/3L66JjlVZDIZZ/UGXInAB4smwzAIDAxETEwMYmNjOdGk/eTTAxJJAeHYWlhc\nXAylUok5c+a4fFufECLJsfD390d0dDSX3DKbzejo6EBdXR3Ky8shFosRFBSE6OhozJgxAxEREYiM\njCTRFCgkkgJBp9MhJycHDQ0No24tJCZOYGAgAgMDR4imTqdDWVkZJ5qOSJNEU1iQSPIcq9WK/Px8\nnDlzBizLIiMjw21/nEKdbk8UZ6LZ3t6OmpoaXLp0yaWu7YT3oU+Nx3R0dCAnJwelpaVcOYs7IZF0\njkM0Y2JiAMClru2E96FPiacM3lo4Z84cj+wmYRgGLMu6/Tp8Z7hr+2DRLCws5EQzPj4e0dHRCA8P\nR2RkJNnC+ShklcYzLBYLzp07h7y8PIhEIqSkpHhsa6Fer0ddXR21UpgiV3NtnzFjxpitLoSKyWTC\niy++iM8++wy9vb2IiorCU089hZtuugkA8Omnn+K5555DZ2cnVCoVHnvsMdxzzz2uuDRZpQkJx9ZC\nrVbLrXd5mumQ3XY3crmc6/8z2LX90qVLuHDhglPXdo1GI1jRrKurw+bNm3HfffdxyxM7d+7ELbfc\ngg8++AC9vb04cOAAcnJyEBwcjPvvvx/btm1DZGQkrr/+erePjyJJnlBSUoKTJ0/CZDIhJSXFK7tB\n6uvrUVtbi5UrV3r82tOFqxkQKxQKJCQkCM61fWBgAMuWLcOzzz6LzZs3c8fNZjPkcjk0Gg2Cg4Nx\n4cIFKBQKPPbYY/j9738PhmHw6quv4tFHH53qECiS5DuOrYUXLlxAQEAAFixY4LU/Dr5uS+QTzgyI\nHaJ58eJFnD9/3mmrC41G4+WRT4433ngDWVlZQwQSALeE1NraiocffphbczcYDGAYBnK5nJuKA1eW\nMObNm4c///nPSElJQWhoKEJCQlzyt0KRpA/T1NTENcZKSEhw2dbCydLY2IiqqiqsWrXKq+OYzgx2\nbXdEmoGBgQgJCUFCQgLvXNtnzZqFEydOjKjMyM/Px5IlS8AwDIqKijB37lwAQGdnJ77++mssX758\nyGt++OEHp21Fhq/Xnz59GllZWYMPUSTJR1iWRWFhIb799luYzeYpOfe4EioB8j5Ccm23Wq3YtWuX\n09K13NxcAEBYWBgnkAAQEhKC2267bcTztVotgJHf0cHVGJs2bRoukOOCRNLHMJvNOHHiBAoKCty2\ntZAQDs5c2x2ief78eU401Wo1Z3TiK6IpkUhw1113OX0sJycHALB69epxnUur1eK9997DnXfeCZlM\nNuRvpq+vD9dddx3+8pe/TG6ck3oV4Rbq6upw/Phx1NfXIykpyWenTHa7naJKH2W0Vhe5ubm8aHVh\nt9tx8uRJMAyDNWvWjOs1VVVVePDBB52auTz55JP4xS9+wRmYTBQSSR+AZVmua6HVasWCBQt8ct8v\nCSP/GM213SGavubaXlBQgPb29gmJ5J133uk0qDhx4gR0Oh3+8Ic/THo8JJJexmQycV0Lw8PDERcX\n5+0hjQnLsrzsjUNc3bW9vb0der2eMyB2fBe94drumGprNJoh65EOtFotjEYjli5dyh1zVi9psVjw\n6KOP4uuvv57SeEgkvUhVVRWOHz+Oy5cve2xr4VSgSFJ4XM21vbW1FbW1tQDc2+riyJEjCA0NxZIl\nS7hjDpG82s6up59+Gps2bRoiks547bXXsH79eu7eJguJpBewWq04d+4cvvvuOzAMg4yMDEilUm8P\na0wcIkmRpHC5mmt7c3MzampqALjOtX379u148803IRaLUVpaiqSkJHR1deHYsWNgGAZpaWkjXqPV\napGbm4u9e/eOeu7W1la89NJLqKysnPC4hkMi6WGMRiOOHTuGioqKSXct9BYUSU4/JuPaHhERMS7R\n3Lt3LxiGgZ+fHzeLev3116HRaFBfXw+r1Trk+f39/bj77rvxyiuvICAgYNRz//Wvf3VZ2xISSQ9S\nWlqKEydOoKOjA3PnzkVQUJC3hzQpaNfN9GW4a3tvby/XureiomJCrS5SUlJQUVGBjz/+GOHh4di3\nbx8OHTqEvLw8rF27Fh999BEeeughxMTEoKysDD//+c+RnZ2NO+64Y9QxsiyLd955Z8QunslCIukB\n+vr6cObMGeTn53NbC/k4XaVIkhhOQEAAAgICxmx1Mdi1PSoqClKpFAcPHsTjjz+Ohx56CAMDA1i1\nahWOHDkCjUaDU6dOYceOHVi+fDkkEgliY2Px5JNP4sYbbxxzTOfOnUNdXR1mzZrlknskkXQzTU1N\nOH78OKqrq5GQkDDlRWRvMnhNkiCcMZpre3FxMbq6uhAfH49f/vKXiImJwUcffeT0PCqVCu+8886k\nxuBY03TVNl4SSTfh2Fp4+vRpmM1mzJs3zye2FhKEJ3GIplqtRkFBAXp7e93emvfbb78FcGULoysg\nkXQDg7cWhoSEID09XVBbC2lNkpgI3d3dKCgogN1uR0xMjNtL3erq6gDAZUEJiaSLqaurQ05ODvR6\nvU9vLZwMtCZJTBSj0Yji4mLIZDKkp6dDp9O5fTfZ5cuXAcBlPYRIJF0Ey7LIz89Hbm6uT28tnAok\nksREMBgMKC0tRWhoKObPnw+pVAqbzeb2mmDHTMdVJh4kki7AZDLhxIkTKC4uhkajQUJCgreHRBBe\npaamBjU1NQgPD8fcuXO55SabzQY/Pz+3Xvuf//wnamtrkZGR4ZLzkUhOkZqaGhw7dgwtLS2YPXu2\nyxaLfRHKbhNjwbIsysvLYTAYEBcXh6SkpBHPcXf526pVq1xqDE0iOUmsViu+++47nDt3jldbC6cC\nTbeJ0bBarSguLkZ7eztmz57t1EzXbrfzrt84v0brI/B5a+FUoEiSuBoWiwUXLlxAX18f5s2bd9Ua\nRbvdzruNFCSSE6S0tBTffvst2traeL21cDJQJEk4o6enBxcvXgTLssjIyBg1YUKRpICxWCw4ffo0\n8vPzIZPJkJmZybtfRFdBdZKEg46ODhQWFnJ/E2PVJlIkKVBaWlpw9OhR1NTUIC4uDhEREd4eEkF4\nnebmZpSWliI4OBjp6enjKnmjSFKAFBYW4tSpU+ju7vaZroXehiJJQqfToaqqiutmON7okCJJAWE2\nm3Hy5EkUFBRAoVAgIyNDUFsLJwOtSRIAUFZWBoPBgNjYWMyZM2fCr+dbJCnov3qTyYSnnnoKycnJ\niIuLQ1ZWFj7//HPu8U8//RTp6emIj49HZmYmPvzwQwBAfX09PvnkE+Tn5yM+Ph5JSUnTXiABym5P\nd6xWKwoKCmAwGDBr1qwJC6TNZgPDMLwTSX6NdgLU1dVh8+bNuO+++1BUVASZTIadO3filltuwQcf\nfIDe3l4cOHAAOTk5CA4Oxv33349t27ahs7MTVqsVFosFGRkZbt8dwCcokpy+WCwWFBQUwGw2Iy0t\nbVKWf3wVSWYC60u8WYgaGBjAsmXL8Oyzzw5xJzabzZDL5dBoNAgODsaFCxegUCjw2GOP4fe//z0Y\nhsHGjRuRnZ1NWwud0Nvbi++//x5paWmCMu4gRqenpweFhYWwWq2YP3/+pPdEm81mlJeX4+6770Zs\nbKyLRzlpxvzlF+Qc8o033kBWVtYI+3bHgnFrayvuv/9+zrLJ0SxIJpNh/fr1VxXI5uZm3Hvvve4b\nOEH4GB0dHcjPzwcALFq0aEqmEXa7HSKRiHeRJL9GO07efvttnDhxYsTx4uJi2O12MAyDrVu3wmq1\nIi8vD1lZWQgNDcV1113ndPdMT08Pvv32W+zZswddXV0euAPfxDHdpuz29KClpQWXLl2aUInPaNhs\nNhJJX8BqtWLXrl1O943m5uYCAMLCwhAZGYnPP/8c5eXliI6OHtL3dzA7d+5Eb28vUlNT3TpuPkBr\nktOHuro6aLVaqNVqzJs3zyVlO1arlZdrkvwa7TiQSCS46667nD7maHqemZmJTz75BEajEcnJyZDL\n5Vc937PPPjvk9dM5knRAkaSwqaioQENDA2JiYjB79myXVXbQdNvHsdvtOHnyJBiGgb+/P/r6+njb\ntZAg3AHLsiguLkZbWxsSExNdnrzka3abX6OdAjk5OWhvbwcALF26lKbPk4DWJIWLxWJBYWEhuru7\nkZqaisjISJdfg6bbPkxxcTHeeustAIBCocCiRYtGPKehoQFdXV1ISUnx9PB4A61JChOz2YyCggJY\nrVakp6e7rbzLsSWRbyIpyBKgI0eOIC8vD319fTh8+DD+/e9/c2U+8+bNc/qavXv3Qq/Xe3KYvIV2\n3AiHzs5OrsRn4cKFbq1/tdlsvFzeEpxIbt++HevXr8eKFSvw1ltvIT8/H2FhYaioqADDMIiPjx/x\nmoaGBly6dAmrV6/2woj5A0WSwqKlpQUXLlxAQEAAMjMz3e6NarPZeBdFAgIUyb1794JhGIhEIphM\nJqSnp+PEiRNc7xmbzTbk+RaLBS+99BJ++tOf0hbEMaC928JBr9ejpKQEoaGhyMzMhL+/v9uv6YlO\nie5AUCLZ3d2NyMhIyGQyPPjgg7jmmmuQm5uL3NxcvPnmm4iNjUVOTg5aW1sBXKkFe/rpp7Fs2TKs\nWbPGu4PnARRJCgOtVgutVovIyEikp6d7bArMsiwv2yzzL/a9CjqdDkePHsXNN9+M06dP4+DBg/j4\n448xb948/O53v0NISAhee+017NmzB4888gjEYjE0Gg1uv/12LF261NvD5wUkkvyGZVmUlJSgtbUV\nCQkJSExM9Oj1+Trd5t+InaDT6fDWW2+hv78fSUlJeOSRR5y2dlUoFNi+fbsXRkgQ3mVwiU9ycrJX\nmtdRJOlFGIaBn58flEolmpqa0NjYCLFYDKVSCZVKhbCwMF5+OL4G1UnyE7PZjMLCQlgsFreW+IyF\n3W7n5ZqkIERSIpFALpdjzpw58Pf3R0dHB4xGI4xGIyoqKlBeXo7g4GCoVCqoVCqEhoaSiS4xLTCZ\nTCgoKIBIJMKiRYu82t3TZrPxMjkqGJFkGIZzGXGIIQD09fVxgmkwGKDX68EwDBdlajQaBAQEjHmN\nnp4e9Pb2AgAuX7581b7C0wHKbvODy5cvo6SkBIGBgUhPT/dIBns0+NjfBhCQSIpEohHlPQDg7++P\n6OhoREdHg2VZmEwmTjS1Wi0qKioQFBSE0NBQhIWFQalUDoky//CHP6CiogL19fUYGBgAAPzP//wP\nYmNjER4ejueee85j9+kLUPKGHzQ0NKC8vBwqlQrz58/3GXGixI2XcESSY0U4IpEIoaGhCA0NxcyZ\nM2GxWGA0GtHW1oa2tjY0NjZCJBJBoVBAo9FArVbjkUce8dBd8Adak/RttFot9Ho9oqKikJyc7DNL\nSxRJehGJRAKxWOw0khwNmUyGyMhIbjN/V1cXWltb0d7ejqqqKlRVVUEmk0GlUkGtVkOlUvHyl9CV\n+MofHDESlmVRWlqK5uZmJCYmerzEZyz42HMbEJBIOtYkp4JcLodcLkdiYiIGBgbQ3t6OtrY2tLe3\nw2AwALhSRqRWq6FWq7n2D9MNiiR9D4vFguLiYphMJqSkpHilxGc8UCTpJRyRpCsTClKpFOHh4VxX\nOLPZjNbWVhiNRtTV1aG2thZisZiLMNVqNZUZEV6hr68PFy9e9HqJz1iwLEuRpLcQiUSTmm5PhMDA\nQMTFxSEuLg42mw0dHR1oa2uD0WhES0sL7HY75HI5lEolNBrNlBom+TLjWfslPIfJZEJhYSEYhkFm\nZuaoLvu+AEWSXkQikbhVJAfjiCDVajWAK1Fme3s7Wltb0djYCL1eP6SYXa1We738ghAebW1tKC4u\nhr+/PzIyMnz+O0Zrkl5GKpVydYyeJjAwEIGBgYiJiQHLsiOK2e12O4KCgoZEmXxOgNCapPdpaGhA\nRUUFQkJCMH/+fJ/fycLXTomAgERSJpOhp6fH28MYtZi9ubkZ9fX1XCmSIxoNDAz08qgJPlFdXY3a\n2lpEREQgNTWVFz+4fO1vAwhIJCUSCVfs7UsMLmYHrqwhORJAWq0WlZWVCAgIgFKphFqthlKp9Okv\nEsMwFEl6CUeJT0tLCxISEjBz5kxvD2ncWK1WLnfAN3z3r3GCyGQyXiQUFAoFFAoFZs6ciYGBAa7E\nqK2tDQaDAQzDDCkzCg4O9vaQh8CHqEWIWK1WFBYWorOz02suPlPBbrdTJOltpFIpL0RyMFKpdEQx\nuyNjXlNTM6KYXalU+kSZEUWSnqWvrw8FBQXo6+tDeno6lzDkE7Qm6QP4+fl5LLvtLhzF7AkJCbDZ\nbNxapmM90263Q6FQcMYc3ihmp73bnqW7uxsXL14EAGRmZvJ2AwNf28kCAhJJX8/uTRSxWIywsDDO\nbchsNnN7zPV6PXQ6HVdmpFarodFoPBZl8i1i5ytGoxFFRUXw9/fH/PnzeZ3gs9vtFEl6G1fvuPE1\nHGVGM2bM4MqMHPvMW1tbUVZWhqCgIE4wFQqFW9YPKZL0DA4Xn9DQUF6U+IwFZbd9AD6++ZNleJlR\nb28vNy13FLM7yowcbka+XmhM/IeamhrU1tYiPDwcc+fOFUSyjKbbPoBYLJ62CYWAgADExMRwxeyd\nnZ1c1txRzB4YGDgkATTZPzyKJN0Hy7IoKytDU1MT4uLikJSU5O0huQyHTRqJpBeRSCTTViQHIxKJ\noFQqoVQqAVxxh3FMy4cXsztEc6KW/kJe1vAWVqsVxcXF6OjoQHJyMmJiYrw9JJfC106JAImk4JHJ\nZCOK2R0JoOrqalRVVcHf35/bATSWZyZFkq7HYrHgwoUL6OvrQ1pamiBbg9hsNt6uqwpGJKfzdHsi\nOIrZh3tmOnoAMQyDkJAQKJVKhIWFOS1mp/fZdXR3d6OgoAB2u53XJT5jYbPZeJudF4xIDm7hIISF\nbk8w3DOzp6cHbW1taG1thU6nQ01NDfz8/LgEkFKppEjShRiNRhQXF0MmkyE9PZ23IjIeWJalSNLb\nDG4GRiI5OYKCghAUFDTEM9OxntnS0gIAaG5uhkajQVRUFEJCQrw8Yv5iMBi4VscZGRm8FZDxQtNt\nH0AsFnMtHPj6YfgSzjwzHWuZzc3N+OGHH4Z4ZoaFhfnElkk+oNPpUFVVJagSn7GgSNIHcEcLB+I/\nOIrZW1tbIZVKER0dPcQzs7y8nCtmV6lUvPfMdAcsy6K8vFyQJT5jQSLpA4zWe5twHWKxeFTPTIPB\nAL1eD4ZhhjizC3W9raenB3//+99x5swZ9Pf3Q6VS4Y477sCKFSsAAKdOncK+ffvQ3d0NiUSCJUuW\nIDs7GzNmzPDyyD2LzWbj7UxDUCLpio6JxNgMz24P9sxkWRYmk4kTTa1Wi4qKCgQFBQ0xGhZClNnS\n0oJnnnkGP/rRj/DnP/8ZUqkUH3zwAXbt2oXHH38c/f39OHHiBJ5//nmUl5fj888/x/Hjx7Fq1app\nJ5IAP/vbAAITSZFIRNNtNzNWdttRqB4aGoqZM2fCYrHAaDRy65nDPTM1Gs2Ei9l9AavVimeffRb3\n3Xcfli1bxh2/7bbbsH//frzzzjsICAjAK6+8gvLychw/fhxVVVUQiUTQ6XRIT0/34ug9D187JQIC\nE0mKJD3DRH6IZDLZCM9MR8a8pqYG1dXVkMlkXP+fsYrZfYVDhw4hJSVliEAC/zElNplMWL9+PUpL\nSyGVShEQEACRSISAgAAsXbp0yGs+/PBDKJVKLFiwAGq1GizLorm5GXl5eQgLC8PatWs9dl/uhCJJ\nL+MQSavV6u2hCJqptm9weGYOL2bv6OhAU1MTAAxxZvfV4up//etfePnll0cc1+l03Psjl8sRGBiI\n9PR0ZGSH2u8YAAAgAElEQVRkYMWKFZgzZw40Gs2Q11RUVODSpUsjzpWamorf/OY37rkBD8PXTomA\nwESSstuewVU7boYXs5vNZq7/T11dHWpra7lSJEcCyBcW/202G7Zt28aNezCXLl3ierCnpKRg7ty5\nEIvFkMlkXDLHGcHBwVwju+joaKxfvx4bN24UTPG+w+CCjwhGJB1Nhmi67V7c+UcbGBiIuLi4IcXs\nji2TLS0tXGted3tmjoVYLMZ1113n9LHTp09DKpUiLS0N8+fPH9f5GIbB3r170d3dDZlM5hM/BK6G\nIkkfQSKRkEh6AE/s3XZWzO6Ymjs8MwcXs3vbM9Ph4lNbWws/Pz8sX758wufwtaZvroQiSR/Bz88P\nXV1d3h6GoPHW9M9RzO7wzOzo6BhSzO6IMh0JIE8Ws1ssFhQUFKC2tpbzDkhLS/PItfkAnzslAgIT\nST52TOQj3nYBGl7M7vDMdDRMG+6ZqdFo3FbM3tPTg8LCQlitVk4M5HK50zpIg8GA7u5uzJ49e8hx\nu92O5uZm/O1vf4PJZMLAwACio6Nx9913c/fIZxw/HCSSPoBUKqXptpvxxUSCM89Mh2hWVVVBq9Ui\nICBgiJvRVP9g8/PzAVxZBpBIJFi0aBF2794N4EpW2hn79+/HwoULR4gkALz33nv46U9/ynlJ7t69\nG08//TReeOEFzkCZr/C5vw0gQJGkSNL9eDuSHAuHZ+bMmTMxMDDAtbIwGo1oamoaUsyuVqsnvBb4\npz/9CYcOHYJEIsH27duxZs0a2Gw2FBUVgWEYxMbGjniNwWBAWVkZHnrooRGPLVmyBNdffz0CAgK4\nY/fffz9+8pOfYO/evXj00Ucn/ib4EFarlUus8hESSWJC8G07oVQqHVHM7siY19TUoKqqCjKZbEj/\nn7Gyy1999RUkEgmCgoKwaNEiyGQyfPzxx1AoFGhraxsxmxkYGMDu3buxbds2p+fesmXLiGNyuRwx\nMTE4e/Ys7rvvPl5Hk7Qm6UP4+fnRdNsD+HokORqOYvaEhATYbDYu+eNYz7Tb7VAoFJxoDvfMrKio\ngEajQW9vLx5//HGo1WqcOnUK58+fx0svvYSdO3fi9OnTuPHGG6FSqdDQ0IB33nkHS5YswcqVKyc0\n1uDgYNjtdhQXF2PVqlWufBs8isPjlUTSB+Drh8AnfHFNcrKIxWKEhYVx64CDPTP1ej10Oh1XZqRU\nKrntlA8//DCOHz+O999/HzabDampqdi5cyfkcjl27dqFffv2YceOHVwZ09atW7Fw4UKnYzh8+DDO\nnDmD7du3X7W3jdFodNt74An43E4WEJhI0o4bzyDU99hRZjRjxgyuzMghjDqdDv39/bjmmmsQGRmJ\njIwMp+eQy+V48MEHx33Nr776Co2NjTh79iyys7OHPGY2m7lx8Rk+t5MFBCaSfP0Q+ISQIsnRGF5m\ndOnSJbS1tXFrm65CoVCgq6vL6e6cjo4OAOC9OS9lt30I6phIuAuRSOQWZ+2MjAzcf//9SExMHHK8\noaEBHR0dSEpKGvEY33Bkt/kqkvxKVY6BwwmIcB/T9f212+1u+QFev3499u/fj/7+/iHHjx8/DplM\nhp///Ocuv6anoey2DyGRSAS7XuZLTNdo3R11fnK5HLfddhteeOEFXHvttVAqlSgqKsK5c+fw7LPP\nIj4+3uXX9DR8b84nKJGk6bb7mc6RpLtISkrCU089hfPnz0Ov1yM1NRX33HMP72pSrwaJpA/hmG47\n9ooS7oF+iFyPn58fr2shR8Nms/E6Qy8oJaEWDu5nuv742O32aXvvU4XP7WQBgYmko90ptXBwL9Nx\n3dddiZvpAN+n24ISSZlMRh0T3cx0XpOcrvc+VSiS9CEcbWVpuu1eKKIiJgKJpA9BHRMJd0JrkpPD\nZrPBz8/P28OYNIL61B2RJEU6hKuhJZypwddCckCAIknZbfdC0RQxUViW5a3hLiBQkaTptnuZjlEV\nlQBNDYokfQSJREK7btzMdM3wTscfBlfB53aygMB23ADUe9sTOPsR6u3txcGDB3Hu3DlYLBYolUrc\neuutWLJkCQDg7Nmz+OSTT2A2mxEcHIzs7GysXr3a00OfEtP1B2Kq2O12XkeS/B35VZBKpRgYGPD2\nMASLM6G4fPkyfve732Ht2rV4/fXXIZVKceDAAbz88st4+OGH0d/fjzNnzuDZZ5+Fv78/3nrrLbzx\nxhsIDQ1Fenq6F+5ictB0e3LwXSQF96lTW1n3MziStFqtePnll3HXXXdhy5YtXD3c1q1bAQB79uzB\nwYMH8cQTT0Aul2Pfvn04ffo0GIaBTqfzyvgJz0LTbR9DJpPR+pEbGR5Nffnll5gzZw4WL17s9Hld\nXV3YtGkT1y61vb0dDMMgICAAS5cuHXH+nJwcnDp1Cn19fTCZTFi0aBFuueWWEQ25PA3LsjTdngSO\n943PUbjgRJIiSc/y1Vdf4bnnnhtxvK6ujtvKN1gMf/aznyErKwtz5syBRqMZ8poDBw7AbrfjmWee\nAQDU19fjueeew6lTp/Db3/4WcXFx7r2ZUaBk4ORwOHLRdNuHoN7b7mVwNGWz2XDnnXeOEDsAKC8v\nB3Clh8uMGTO444GBgVixYsWI12i1WlRUVODOO+/kjsXGxuJnP/sZTCYTXn31VVffyoShSHLi8L2/\nDSBAkaTptmdwrDNdzQOxpKQEADB37txxne/IkSPYsGHDiOMLFy5EUFAQDAYDioqKJj/gKUKR5ORw\n9Nzm85qk4ERSKpXSF9qNjCeastvtKC4uBsMwSEtLG9d5q6ursXv3bly4cGHE9aKjowFcmcJ7Ez6v\nq3kLvve3AQQoktTnxjOM9kNUW1uLnp4eABi3SNpsNlgsFuTn5zt9bKxruhv64Z0cjkiSRNKHEIvF\nJJJuxBFJjiYaxcXFAK40uRq8HunAYDCgsrJyyLH/+q//wuzZs7Fu3bohx202GxoaGgAAM2fOnNLY\npwKJ5OSwWq28jyT5O/KrwOcPg084RKOgoABBQUFISkriHnOIZGpqqtPX7t+/HwsXLsTs2bO5YytW\nrMCKFStGPPfMmTPo7+9HXFzcuNc33QHt3Z4cjrVrPv9d8nfkV4H2brsXh1DY7Xa8//77+PLLLyES\nibB7925ERUWht7cXRUVFYBgGsbGxI15vMBhQVlaGhx56aMxr2e12HDp0CGKxGA888IDL74VwP5Td\n9kFkMhmJpIfIyckBwzCQSqVcN7wvvvgCCoUCAEbUqw4MDGD37t3Ytm0bZDLZmOf/4osvoNfr8eCD\nDyIlJcX1NzBBqARo4litVlqT9DUoknQvg9ckY2JiEBAQgCeeeAIhISE4deoUzp8/j5deegkxMTE4\nffo0jEYjAKChoQG7du3CkiVLsHLlyjGvU1NTg7///e944IEHsHbtWrfe03ig6fbkEEJ2m78jvwp8\n/jD4gtVqBcuyeOKJJ/B///d/ePfdd2Gz2ZCamoqdO3dCLpdj165d2LdvH3bs2AGxWAy1Wo2tW7di\n4cKFY56/s7MTr7zyCn7yk5/4hEAClLiZLHzvlAgIVCQd7uR8LmD1RViWRV1dHdra2pCXlwe1Wo07\n7rgDarV6xPRZLpfjwQcfnPA1rFYrXn31Vdx+++249tprueO9vb3o7u5GWFjYlO9jMpBITg6bzTau\npRVfhkSSGBf9/f0oKSmB2WxGVlYWbDYbjEYjtFottFotgoODoVQqoVarIZfLJ32dv/71r1i9evUQ\ngQSuZNH7+vqwZs2aKd7J5CCRnBwUSfogYrEYIpEIVquV979gvkJ3dzcuXboEm82G9PR0LjGTkJCA\n/v5+GI1GGI1GNDY2Qq/XQyqVcoKpVCrH/WP15ZdfIjo6GuvWrRtS68qyLIqKirwmkA5oTXLi8L2d\nLCBAkaSOia7FaDSivLwcUqkU6enp8Pf3H/K4n58foqKiEBUVBbvdjvb2dhiNRrS3t6OlpQUikQhy\nuRwqlQpqtZqzTBtOYWEh9u7dC5Zl8be//W3E4wzD4J577nHLPY4Xym5PHIokfRCHSJJd2tRpbGxE\ndXU1FAoF0tLSxowIGYaBSqWCSqUCAJjNZrS1tcFoNEKn06Gmpgb+/v5Qq9VQq9VQKBSc8Lz77rtc\nJtQZoaGhVxVYT0A/upODIkkfhDomTh273Y7q6mo0NjYiMjISSUlJk4qiAgMDERgYiBkzZsBqtcJo\nNKKtrQ3Nzc1oaGiARCJBaGgoVCoVXnvtNZ/+YyLT3cnBsizvl70EKZI03Z48LMuitLQURqMRCQkJ\nTvdeTwaJRILw8HCEh4fDbrfDZDJxa5mtra0AMGRaHhQU5JLrEt7Fbrf79I/feBCsSNJ0e+IMzmCn\npqY6NdN1BQzDICQkBCEhIUhMTER/fz/a2trQ1taG+vp66HQ6+Pn5QalUctN3X4jiKHEzcViW5X3t\nMr9H7wSHSNJ0e2JcLYPtCfz8/BAdHY3o6GiwLMslfzo6OtDU1MSJqkajgVKpHJE88gQ0M5kcLMvy\nvhRPkCLJMAx9qSfAWBlsTyISibjEDnBFvB1rmVqtFsCVtU7HtNyTYu4L0SwfoUjSB6FmYOOnsbER\nVVVVCAkJGVcG29MEBwcjODgYcXFxGBgY4LLlTU1NqK+vh0Qi4WoyVSqV28ZPe7cnB9/byQICFkmL\nxeLtYfg0rspgexKpVIrIyEhERkbCbrejo6ODS/5cvnwZDMNAoVBAqVRCo9G4tGTI198bX8Vut1Mk\n6YtIpVL09vZ6exg+i7sy2J6EYRgolUoolUrMmjULvb29aGtrQ3t7O3Q6HWpra+Hv74/Q0FBu589k\nhc6x+4eEcuKQSPooNN2+Op7KYHuagIAAxMbGIjY2FlardUiU2dTUBJFIxNVkOjPkGAta454cNN32\nUaj3tnO8mcH2JBKJBBqNhvsB6Orq4tYyJ2PIQd+lyeEowKdI0gehSHIkvpTB9jRyuRxyuXzShhwO\nkaTEzcQQQqdEQKAiKZPJSCQH4esZbE/iKkMOYmwc/W34/n0TpEj6+fnRGhL4mcH2JFcz5HAkfxyG\nHMHBwTCbzV4eLf9gWZYiSV+F7x+KKxBCBtvTDDfkaG9vR1tbG1paWqDX6yEWi5GQkMAJK9/3JLsb\nIXRKBAQqkmKxeFovtgs1g+1JJBIJwsLCEBYWhvj4eJhMJoSHh8NsNpMhxzihNUkfhu8fylSYLhls\nT2K32xEUFIS4uDiEhYXxxpDD21B224eZrm1ljUYjysrKIJPJpl0G290M3pY4XkMOx1bJ6fo5OPpM\nkUj6IBKJZNqJJGWw3Yfju+QsOhzNkKOqqgpVVVUICAgY4sY+XXCsSfK9dEqQIjmdIknKYHuG8b6n\nww05HIJpMBhGGHKEhoYKOvlD2W0fhu8fynihDLZnmOwPrlQqRUREBCIiIjxuyOEL0JqkDzMdem9T\nBttzuGJW4klDDl9BCJ0SAYDfiwVXQegtHLq7u1FQUID+/n6kp6eTQLoZh0i6cm3NYcgxf/58LFu2\nDKmpqQgJCYHRaERJSQlyc3NRUlICg8HAW9s/Z50STSYTnnrqKSQnJyMuLg5ZWVn4/PPPucc//fRT\npKenIz4+HpmZmfjwww89PewRCDKSFIvFXAsHvndqG05bWxvKy8spg+1B3CGSg3G1IYevMDySrKur\nw+bNm3HfffehqKgIMpkMO3fuxC233IIPPvgAvb29OHDgAHJychAcHIz7778f27ZtQ2RkJK6//nqv\n3YcgRVKoLRwogz09uJohR0NDw7gMOXyFwZHkwMAAtm7dihdeeAGbN2/mnvPrX/8azz//PH71q18h\nODgYFy5cgEKhwGOPPYb9+/eDYRgUFRWRSLoaoU23KYPtXUYrAXI3fDbkGBxJvvHGG8jKyhoikAA4\ngW9tbcXDDz/MlUgZDAYwDAO5XI6bbrrJ6fl/+OEH7Nq1CwaDgYvG33zzTZcnMAUrkgzDCKJjImWw\nfQNfmJWM15DD8ZzQ0FCv/pgOjiTffvttnDhxYsRziouLYbfbwTAMtm7dyh3/05/+hJtvvhnLly93\n+p3//PPPcd9992HPnj24+eabAQBff/01br/9duTm5rr0PgQrkkLYv00ZbN/AV/0knRlyGI1GtLS0\noLGxERKJBCEhIVyU6elMM8uykMlksFqt2LVrl1OxcwhaWFgY5s6dyx0PCQnBbbfd5vS8JSUluPfe\ne/Hqq69yAjkwMIAnn3wShYWFMBqN3A+JKxCsSDpKgPjK4D3YGRkZvFmsFyKOSMeXlzgGG3LY7XZ0\nd3ejtbUVRqMRlZWVqKys9Lghh91uh1QqhUQiwV133eX0OTk5OQCA1atXj/u8//u//4vo6Gg88MAD\n3LGBgQE0NDQgJibG5X8rghRJmUzG6zVJymD7Hr4w3R4vjrU8uVyOxMRErxlyjNUEzG634+TJk2AY\nBmvWrBnXOQsLC/H1119jx44dQ44HBgaiuroaDMO4PGIWpEiKRCLeRpKUwfY9+L5s4y1DjrE2cxQU\nFKC9vX1CInnw4EEwDIP58+ePeCw4OHiyQx0VQYokwL8+N5TB9m2E8lmM15DDMS0PCQmZ0vVGiyQd\nU22NRjNkPdKBVquF0WjE0qVLR7wmIiICn3zyCd577z1cvnwZarUazz///JDnugrBiqRMJkNfX5+3\nhzEuKIPt27i7mNybODPkcLThbWhomJIhx+B2skeOHEFoaCiWLFnCPe4QvGuuucbp659++mls2rRp\niPBVV1cDAE6fPg2FQoGvvvoKAPDWW29h1apVOHz4MK677rqJvQljIGiR5ENfEspg+z58Wo+cCq42\n5HCsSW7fvh1vvvkmxGIxSktLkZSUhK6uLhw7dgwMwyAtLW3Ea7VaLXJzc7F3794hxy9fvgwAOHHi\nBI4dO8Yd/+Uvf4k//vGPuO+++1BRUeHSdXzh/TT+//Bhuj14D3ZGRgYJpI8jxEjyajgMOWbNmoUl\nS5Zg8eLFSEhIAMMw0Ol0+P7775GXl4fKykoYjUanPyQOkdy7dy8YhoGfnx9XLP76669z3/fh9cz9\n/f24++678corr4wQYj8/PwDADTfcMOJ6K1asQENDAz799FOXvAcOBPupS6VSn15wb2trQ2FhIUQi\nETIzM6nEx4dxfI+Esi45GQYbcqxYsWJchhyO6XZKSgoUCgX+8Y9/IDw8HPv27cOhQ4eQl5eHlJQU\nfPTRR2hoaAAAlJWVYcOGDcjOzsYdd9wxYhxKpRIAnEafjvInxzTeVQh2uu3LkSRlsAk+IxKJnBpy\ntLe3c4YcgYGB6OnpgVgsxsGDB/H444/joYcewsDAAFatWoUjR45Ao9Hg1KlT2LFjB5YvXw6JRILY\n2Fg8+eSTuPHGG51eW6PRQK/XO3V4d0SZjY2NLr1fwYqkv7+/z4mk3W5HVVUVDAYDZbB5hK/uuPEV\nnBlylJWVcdsSY2Ji8NFHHzl9rUqlwjvvvDPua6WlpeHChQswmUwjHnNM+aeakR+OYD91X2vhwLIs\nLl26BIPBgMTERMyePZsEkhAcfn5+6OnpQVhYGLZt24bExESXnn/NmjWw2+1obm4e8ZjJZALDMMjI\nyHDpNQUrkr5kGd/f34+LFy+ivb0dqampiI2N9faQiElAP2pjo9Vq0dfXhxtuuMHpuuFUufnmmxEY\nGIji4uIRj/3www8Qi8W4/fbbXXpNwYqkrxhcUAab/9B0e3xUV1ejt7cX119/PebNm+eWayiVSjzx\nxBP4+9//jp6eHu54XV0dzp8/j1/96leYOXOmS6/pO+GWi/GFSJL2YBPThZqaGvT09OBHP/oR0tPT\n3Xqt//f//h9qamqwYcMG/OpXv0Jvby9eeOEF/OxnP8OLL77o8ut5X0nchLczxpTBFg4OFyDCOTU1\nNeju7sa6devcLpDAlWWPPXv24Pz58zh37hz8/Pzw2WefITk52S3XE6xISiQSr2S3KYNNTCdqa2th\nMpmwbt06LFiwwKPXXrp0qVv2ag9HsCLpjcht8B7sxMREStAIBIoknVNXV4fOzk6sXbsWCxcu9PZw\n3IZgRdKxJumpLzjtwRYuvpAA9DX0ej3a29uxdu1aLF682NvDcSuCFkmRSASWZd0eVZKLuPChSPI/\n6PV6tLW1TQuBBAQuko5mYO4UScpgCx9f2pTgbRoaGtDa2oo1a9YMsT0TMoIVSalU6vYWDpTBnj5Q\nJHnl+97c3Ixrr70WWVlZ3h6OxxCsSLqzGRhlsKcXLMtO+0LyxsZGNDU14dprr8Xy5cu9PRyPImiR\ndEckSRlsYrrR2NgIg8GA1atXTzuBBAQukgzDuDQzSRlsYrphMBhgMBiwatWqq7ZZEDqCFkmxWOyy\nSNKRwWZZljLY04zpWgLU1NQEg8GAa665BitWrPD2cLyGoEXSVWuSlMEmptuaZEtLCxobG7Fs2TKs\nXLly2t3/YAQtko46yeH09vbi4MGDOHfuHCwWC5RKJW699VaupOHs2bP45JNPYDabIZPJMH/+fCxb\ntgxz586lDPY0hGXZaZWYa2lpQX19PZYuXYrVq1dPa4EEBGyVJhKJnE63L1++jKeffhoKhQKvv/46\n3nnnHSxYsAAvv/wyTp06haNHj+Kbb77Bb3/7WzzyyCMIDg7G119/PWajdULYTJdaycuXL6Ourg5L\nly7FmjVrpr1AAgIWSWBknxur1YqXX34Zd911F7Zs2cL1EN66dSsAYM+ePTh48CAee+wx6PV6fPbZ\nZ6ioqADDMNDr9V65B8L72O32aSEWbW1t0Ol0yMrKIoEchGCn28AVkezt7eX+/+WXX2LOnDkjtlI5\nvgxdXV1Yv349KioqYDabIRaLwTAMAgICnLqNnD17FocOHYJMJoPVakV8fDzuvvtuSuoIjOlgcNHW\n1oaamhosWbIEa9euJYEchKBFUiaTDXEv/uqrr/Dcc8+NeF5dXR33hyCXy2GxWJCRkYGFCxfi4sWL\nmDNnzohyn1OnTuGPf/wjfvvb3yI1NRV2ux2vvvoqnnnmGbzwwgtjNm4nCF/BaDRyArlu3ToSyGEI\n+t0Y7Clps9lw5513Oq1tLC8vB3Clt3B4eDgWLFgAuVyOwMBArFixYsRrenp6sGfPHixbtgypqakA\nrmxb++///m/o9XqXN0cnvIuQ1yONRiOqq6uRmZlJEeRVEPQ7IpPJuOy2WCzGqlWrnD7v+++/BwAk\nJCQgMzNzzBKfM2fOoKura8Tug/DwcMycORNHjx71uXa2xNQQ4nTbaDSiqqoKmZmZ+NGPfuQTLU98\nEUGL5PDEzXDsdjsqKyu5SHL58uXjymDn5+cDgNOGQ4mJiTCbzdw5Cf4jxBIgo9EIrVZLAjkOBC2S\nfn5+VxVJlmVRXFyMoqIi9Pf3g2GYcXd40+l0AK40Vh+OY2peW1s7uUEThJvp6OhAVVUVFi5ciBtu\nuIEEcgwELZJisdjpepKjD3ZnZycsFgsAQC6XY8aMGSOeazAYUFlZyf2fZVkYjUZIpVKnX67g4GDY\n7XZ0dHS48E4IbyOUtbqOjg5UVlYiPT2dBHKcCOOTvwoSiWSISBYUFKCoqAgFBQVcBrumpgYAuATM\ncPbv34+Ghgbu/2azGSzLcjWWw3EcH5xVJ/iNUPZum0wmTiDXr19PAjlOBP0uDY4k33//fXz55Zdg\nGAYPPPAAVq1aBbvdjqKiIjAM49TyzGAwoKysDA899BB3zBF5Xm3t0hFxmM1mV98OQUwak8mEsrIy\nTiBlMpm3h8Qbpk0keezYMe7YokWL4O/vjy+++AIKhQIARqxdDgwMYPfu3di2bduQL9RYC/j9/f0A\nhDM9I/i/46arqwvl5eVIT0/Hhg0bSCAnCH8/+XHgiCS1Wi1CQkLg7++PX//611CpVDh16hTOnz+P\nl156CTExMTh9+jSMRiOAK308du3ahSVLlmDlypVDzjlWeZAj0iSnIMIX6OrqQmlpKebOnUsCOUkE\nPd2WSCRoa2uD2WzGL37xCxw5cgTvvfcebDYbUlNTsXPnTsjlcuzatQv79u3Djh07IBaLoVarsXXr\nVqe9hAMCAhAQEACr1er0mn19fWAYBqGhoe6+PcJDeKLjpjvo6elBaWkp0tLSsHHjRhLISSJokWRZ\nlvOVFIvF2LZtG5RK5Yipk1wux4MPPjju88bGxqKqqsqpM1BnZycAICoqauo3QBCTpKenB5cuXUJq\naio2btxIM5spIGiRXLBgAaKioqDX66HVatHY2Ijq6moEBgYiPDwcKpVqUhHCnDlzuPMNLxtyuAXN\nnTvXJfdA+AZ8WpPs6elBSUkJUlNTsWnTJhLIKSJokQSAiIgIREREYPHixWhtbUVdXR0qKythMBhQ\nW1uLgIAAaDQahIWFjVswly1bhn//+9+4ePHiEJHs7u5GRUUFZs+e7bTQnOAnfNpxYzabUVJSguTk\nZIogXYTgRXIwGo0GGo0GCxcuREdHB2pra1FVVYX6+nro9XoEBARApVIhLCzsqnWQwJWayjlz5uDw\n4cPYsGED99ycnBywLIt7773XU7dEEByDBXLTpk0IDAz09pAEATMBhxPBWqGYTCZOMOvq6tDT0wOZ\nTMZFmM4Es6WlBc8//zzUajU2btyI5uZmHDhwALfccgtuvfVWL9wF4S7y8vIQGBiItLQ0bw/lqjgE\nctasWcjOziaBHD9jThFIJIfR3d2N2tpa1NTUoKamBt3d3ZBKpZxgDs4Q2u125Ofno6GhATKZDBkZ\nGYiOjvbi6Al3cP78eQQFBfmsSJrNZly6dAkzZ87E5s2bERwc7O0h8QkSyalgNpuh0+lQU1OD6upq\nmEwmSCQSqNVqhIWF0XrPNOH8+fOQy+VX3brqTUggpwyJpKvo6+uDTqdDdXU1ampq0NnZCbFYDJVK\nBY1GQ9MbAfPdd99BoVD4nEj29fWhuLgYCQkJyM7OJoGcHGOK5LRK3EwFf39/JCcnIzk5GX19fdDr\n9aitrYVWq0VpaSkYhoFSqURERAQJJuF2BgskRZDuhSLJKWKxWIYIZkdHB+x2O0JDQxEREYGgoCBv\nDwanhDkAAAp5SURBVJGYImfPnoVKpUJycrK3hwLgPwIZFxeH7Oxszn+AmBQ03fYkVqsVer0eer0e\nZWVl6OjogNVq5QSTuijyE18SSYvFgqKiIsyYMQObN2+m7a9Th0TSW1itVjQ0NECv16O8vBxGoxED\nAwMIDQ1FeHg4/frziLNnz0KtVmPOnDleHYfFYkFhYSHi4uJIIF0HiaQvwLIsGhsbodPpUFlZicuX\nL8NisSAkJAQajYZ25/g4ubm50Gg0XhVJRwQZHR2Nm266iQTSdZBI+hosy6KpqQk6nQ5arRYtLS3o\n6+uDQqFAWFgYQkNDebVPeDpw5swZREREICkpySvXHxgYQFFREaKiorBlyxb6UXUtJJK+DMuyaGlp\n4QSzqakJfX19CA4ORnh4uFPHIsLznD59GpGRkV4RyYGBARQWFiIqKgrZ2dkkkK6HRJJPOASzqqoK\njY2N6O3t5RyL1Go1CaaX8FYk6Yggw8PDkZ2dzXXiJFwK1UnyifDwcISHh2PJkiVXdSxyCCYfTWD5\nCsuyHv+BGhgYQHFxMcLCwkggvQyJpI8y2LHIaDSirq6OcyzS6XScxZtarR7VsYjgHzabDSUlJdBo\nNNiyZQsJpJchkeQBKpUKKpUKCxYsGOFYVFdXB39//1Edi4ip4yk/SZvNhsLCQqjVamzZsgXh4eEe\nuS5xdUgkeYZCoUB6ejrS09M5x6Lq6mrU1taisbHxqo5FxOQZ3knTndcpKiqCSqXC5s2bSSB9BBJJ\nHhMcHIx58+Zh3rx5MJvNqKmp4TLlzc3NEIlE5FjkAhiGcXskabPZUFxcjNDQUGzZsgWRkZFuvR4x\nfkgkBYLDFDYtLY2zeHP4YpaUlJBj0RRwd99th0CGhIRg8+bNJJA+BomkAAkMDERqaipSU1Ov6lik\nUqkQHh5OgjkOJlAmN2FYlkVJSQkUCgU2b95Mps0+CNVJTiOGOxa1t7cDADkWjUFOTg6SkpJGdMac\nKizLori4GEFBQdiyZQtiY2Nden5iXFAxOeEch2NRTU0NJ5gsyyIkJIQci4Zx/PhxzJ4926UiSQLp\nM5BIEmPjcCzS6XSoqKhAe3s7ORYNIicnB3PmzEFMTIxLzseyLC5dugR/f39s3rwZcXFxLjkvMSlI\nJImJwbIs12K3oqICra2tnGNReHj4tHOfYVkWJ0+eRHJyskvWCwcL5KZNmxAfH++CURJTgESSmDwO\nize9Xs85FvX29iIkJGTaOBaxLIsTJ04gJSVlyiLpEEg/Pz9s3ryZBNI3IJEkXMNgx6LKyko0Nzej\nt7cXcrlc0I5FVqsV33777ZRFkmVZlJWVQSKRYNOmTUhMTHThKIkpQAYXhGsQiUSIjIxEZGQksrKy\nRjgW1dTUDDHgEIpguqr8p7S0FGKxmASSh5BIEpNiuGORYz95Q0ODoByLHCI5FdEngeQ3JJLElHE4\nFi1evHhUx6KwsDDeCuZkRbK0tBQAsHHjRsyaNcuVQyI8BIkk4VIGOxZ1dHRwgllXV4f6+nrIZDJB\nOBb19vbi4MGDOHfuHCwWC5RKJW699VYsWbIEwJXmYR988AFMJhPCwsKgVqsxe/ZsL4+amAyUuCE8\ngslkQl1dHedY1N3dDZlMxhlw+KpjUV9fH3JzczF//nyEhYUBAC5fvozf/e53WLt2LTZs2ACpVIoD\nBw7gs88+w8MPP4z+/n58/fXXuO2225CdnY0XX3wR+/fvxzfffIPrr7/ey3dEDIMSN4RvoFAoOMei\n7u5u6HQ6VFdXo6amBs3NzZwBhy86Fg12ALJarXj55Zdx1113YfHixdzxrVu34rPPPsOePXsgEonw\ni1/8AjfddBPeffdd7N+/HwzDoKioiESSh5BIEh4nODjYqWNRdXW1zzkWsSwL4D9rkl9++SXmzJkz\nRCAHP97Z2Yk1a9bgpptuQmpqKgwGAxiGgVwux0033TTm9Q4ePIi//OUv+Oqrr1x8J8RkIZEkvMpw\nxyKdToe6ujqfcSyy2+1DyoC++uorPPfccyOeV1dXB5vNBovFgp/+9KdITU0FAPzpT3/CzTffjOXL\nl4+599tsNuPRRx+lBI+PQSJJ+Az+/v5ITk5GcnIyrr322iGORWVlZQAApVKJ8PBwjzoWOUx3bTYb\n7rzzTqc9Z06dOoX+/n6EhIQgOzubOx4SEoLbbrttXNfZtWsX9Ho9iaSPQSJJ+CQymQyzZs3CrFmz\nOMF0OBZVVFSAZVnOgMOdjkWDo0ixWIxVq1aNeE51dTVKS0thtVpx3XXXTeo6FRUVaGlpmfQ4CfdB\nIkn4PBKJBImJiUhMTMTq1auHOBZVVVXBZrNxBhyediyqqalBd3c3ampqIBKJsGbNmkmd56mnnsJb\nb72FPXv2uHaAxJQhkSR4hUQiQXx8POLj47Fy5UrU19dz/clrampgtVq5/eSucCwanrgZjEMgY2Nj\nYTKZwDDMpETy008/xapVqxARETHV4RJugESS4C0ikQhxcXGIi4vDihUrOMei8vJyzrlILpdPybHo\nanXEtbW1MJlMWLduHU6ePAngys6juXPnjniuVquF0WjE0qVLRzxmNpuxd+9e/POf/5zw2AjPQCJJ\nCAKRSITY2FjExsYiKysLTU1N0Ov1qKysRENDA6qqqhAcHDxpxyLH8wsKCtDV1QW5XI61a9di4cKF\n+O1vfwsAuOaaa5y+9umnn8amTZuciuSuXbvw5JNPCsYQRIiQSBKCQyQSITo6GtHR0UMciyorK2Ew\nGFBdXY3g4GBoNJoxHYsckSTDMHj//ffxr3/9C1arFfv378fixYvR1dWFY8eOgWEYpKWljXi9VqtF\nbm4u9u7dO+Kx8vJyXL58GStXrnTZvROuh0SSEDzOHIu0Wi0aGxvH5VjkEMqjR4/CarWCYRguanz9\n9deh0WhQX18Pq9U65HX9/f24++678corryAgIGDEeZ966im8/fbbbrhjwpWQSBLTiuGORY7e5Hq9\n3qljkSNxYzAYIJfL0d/fj4MHDyI8PBz79u3DoUOHkJeXh7Vr1+Kjjz7CQw89hJiYGJSVleHnP/85\nsrOzcccdd4wYxyeffILVq1dTsoYHkMEFQQDo6OjgBLOurg49PT2QyWSQyWSoqqpCWFgY0tPTceDA\nAeTn52NgYACrVq3Ca6+9Bo1GA6PRiB07duDw4cOQSCSIjY3Fk08+iRtvvHHEtXp6evDjH/8Y//rX\nv0ZM9R1lRMePH/fUrU93qH0DQUyUwY5F5eXlaGxsxNatW7F69WqXnP/YsWN4/PHHR5QoDQwM4MyZ\nMwgNDcWCBQsAAIcPH4afn59Lrks4hUSSIKZCV1cXGhoakJKS4vZr6XQ6JCYmUiTpWcYUSao7IIhR\nkMvlHhFIwnchkSQIH8Gxd9toNLqsARkxdUgkCcLL/POf/8T/194dowAIAwEQtPA7+UP+/yorLTeI\nhaAzderlOEIyxtjmnNfjvOcdT7F8n50k8Gd2kgBPiCRAEEmAIJIAQSQBgkgCBJEECCIJEEQSIIgk\nQBBJgCCSAEEkAYJIAgSRBAgiCRBEEiCIJEAQSYCw3zi7/AsC4GtMkgBBJAGCSAIEkQQIIgkQRBIg\niCRAEEmAIJIA4QBEuILoFg1YnAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fef41f411d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from itertools import product, combinations\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "ax.set_aspect(\"equal\")\n",
    "\n",
    "# draw cube\n",
    "r = [0, 1]\n",
    "Y, X, Z = np.meshgrid(r, r, r)\n",
    "x = X.reshape(-1)\n",
    "y = Y.reshape(-1)\n",
    "z = Z.reshape(-1)\n",
    "\n",
    "X, Y = np.meshgrid(r, r)\n",
    "ax.plot_surface(X,Y,1, alpha=0.2, color='grey')\n",
    "ax.plot_surface(X,Y,0, alpha=0.2, color='grey')\n",
    "ax.plot_surface(X,0,Y, alpha=0.2, color='grey')\n",
    "ax.plot_surface(X,1,Y, alpha=0.2, color='grey')\n",
    "ax.plot_surface(1,X,Y, alpha=0.2, color='grey')\n",
    "ax.plot_surface(0,X,Y, alpha=0.2, color='grey')\n",
    "\n",
    "#for s, e in combinations(np.array(list(product(r, r, r))), 2):\n",
    "#    if np.sum(np.abs(s-e)) == r[1]-r[0]:\n",
    "#        ax.plot3D(*zip(s, e), color=\"b\")\n",
    "#ax.scatter3D(X, Y, Z)\n",
    "\n",
    "s = ['$x_0$', '$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$', '$x_6$', '$x_7$']\n",
    "for i in range(8):\n",
    "    ax.text3D(x[i], y[i], z[i], s[i], fontsize=24)\n",
    "ax.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 3],\n",
       "       [6, 7]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(2, 4) + np.arange(0, 8, 4).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "\n",
    "x, y, z = symbols('x y z')\n",
    "p0, p1, p2, p3, p4, p5, p6, p7 = symbols('p_0, p_1, p_2, p_3, p_4, p_5, p_6, p_7')\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4YAAAARBAMAAABz1q9YAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMkSJuyJ2q1Rm\n75mlXrWwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHCklEQVRYCe1Xb4hUVRT/vd35tzur+4ywjGCX\nMYMsRArBSHEzyP47YCoWtFNCJsW6fioDaz8UfRDbqUBYyWlUiqIPTUmEiTqftCBtFYOkoFFCC2R1\nw6Utxemcc+e+e+97b9Rd/RTeD2/fPef8fuf3O3fem1ngxvq/TyBbjHe4Oz4cjV51oQWdDMaCR2+v\nI+EEqSZYztKveeRh/9dMOAkTmAwmLNzZX0fCCVJNsJxVX/PIHevXg3ASJm6cYfgUZH/lSa547dGt\nvtTO+BOZI1j2U6JKW/lQJNe8NbRdcl6piM1cJoQGo8ulyApHChV5cmH+CSnli1XeaBoij2JiBIU6\nu5hYwkYzAV6ZUDsMUQnaMXHZzk2MO5irHLmDUTIKOyupfrkt7EXbOW9D6TTvhDDV/gDWSC6VHcZH\nfCdWDEaXS5EJRwuFHLcMjaheXG/KNUuIPIqJERTq7GLiCHUzAV6ZUJeHqATtmLhs5ybGHcxVjtzB\niIxMZTY6zrWtLSBTOYl0v4eMTFkIN7WU8SO2bAE2Te3GqwxgKwYDXc4pK+wUYjHnmBw+3uFKWaqc\nkwGLIU+dPhGHUYL63vCNIOZSVKnS9jDGEAL9jaQtGYoQacuhIWRMjDYJcxUt1fmxR35o3pmTccZf\n/DCEsUaeWls1Do2Jl6dNczA33cfrDg8XkBo9gG9J7HF0FoElPlK53MyjuVzZ76xg3e2FbBn+lCL+\nwQu53MJcbrbBkAsqj1C5hU/t5Eko8sR6uOWcVCwO+UvYYWO8mSz23h4RNH2gvbshyOl8BqdtjCsC\n2dEgSZIdQh/7/ThCxsRok7DTubdeDsjZoWOFk0DEeEt+s41xRu7PygzHjfyLep0eMjNJVscrOYr0\n8Dx0+ny3tAfJMgX1j6ReHxtb8yArvX7iPMf5Q2FhkCzTvrEMFQVMIXbRVpHjGV2rQulhSSJZpr2N\nmY+lA3EYETSQGDeCBCidT9E7I4SxREwnE45DAZIvbIS3zo8jBGOQLNPF1qbCEhTK9PCzD9LOdWg6\nSzJqfCtmhDDyHIqgtrJHD1N05B8vX6mPySFMzUXXifOYUgG9LreBHkhZinA1vEtdPfgFOIm2YU6I\nNoPR5YKxwlZhTY5JkdMTbRaXq6RmCchr8zCYbwhyMCyoNZ8ZM4KCzjXQcxjCGEK8R+fhOFRCmBBt\nXXSG2qFFyJioNqFSaEBMLOed69B0lqRjQjA00xDGjJwOg5YWZKgq+DqE4TqgpYZDN4+jtYp0Ga8j\nW0M7hxXhHHQcoQdiAShDaVpCaDC6nFMWlV3oyxkKOfYBZa6UxSwqqVkCchrpfl8JcjEsiMTRbLUg\ni8o7hhDGECa+JIzjUIkQwsf5DKOEgolqk7BCK8/+c0MnmneWZNT4hQjGjLxzd8kSZEzAK4T7iI7O\n/sSxxChai5jan7yIhw58P4/jivBvPNmz1Mf9/BHvrXJcCA1Gl3MKJmwXQj2HTJ5Y8MoHeSnV5Sqp\nWQJySlNPEeRiWBBYbCBI6Lgz8D6FXYwhTHl0ho5DAUIIC3yG2qEhFExUm4QVWnnG2zjYvDMnXROs\nNnMJy0IYM/LBWnKuEWRMoJ10ug5FR++qEn+/0G8W75PSOGodR/McF8Lk2NAf/MV0HLht7dkBjguh\nwehyTsGE7UJ1TEI+FYdmSaVcuFwlNUtATu+oGv0QYkEORgQBfYTXgiwqfpe6GEP4HZ+h41CAQtgx\nwGcYJRRMVJuEBU0XZYKH1KwzKOmYEIx3AZ15F2NG3lXBmBFkTGAJdXT7iI7VHD7PhHTE3RLiixDS\nixug70Mi5NczLyG0MRJVFyccFKpjcshNuU6qiIU5LBFLkOxFEDpqslGC5JY70xqs0MXGBIRekc/Q\nTfJeCGeAz1A75DuIlQZGApZxJyydt6MrND5TjkZSschVMGPorNLOUmtGTodxkStDI1cfXgfDVcAc\nvsjv0h1KCe/pWAfoQt9Y9OLi36Xt3Z4cJG6liIXhXbAMFYVMobxLHXKFkHJJagaDSRaSPfTvhYxG\nZ+mvCMIKPG0ESZap6NtwsBrCBITZkZFL4SQjhXDfyNkzMYQKI/x0CVHpsJj4GYN+s86QpC6Xv4K5\nhx8bx6EZeUuVHpuYkeNOwjsY4aNnmlaJf30twAYJmQs900BbLVtGazVdDOIWJojRTTSsIthLySi5\nlYyyLN7zjR+DEUGJu/esjxM0H9viMI0+GXoOoyKEkB5gP44QjLGXRdUIq8jz/HaPklvJKMtKfBWH\nEUHJ/mx3rKBfiSfaZ9dffDapU/TDalWpx25F79NFd3Ggr4/+KRn63eQsjAnSW9FQNcIq8u6iz2LI\nG+WSjLL8Vq/HYJSgKfX6uThBD2/5NAbT6IOD/1aaOux48/M4QsFEtdnhhue+QvPOKU7aS2E6Ljfy\nw/Ejp/el3ec/q4PvJcXUuBwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- p_{0} y z + p_{0} y + p_{0} z - p_{0} + p_{1} y z - p_{1} z + p_{2} y z - p_{2} y - p_{3} y z + p_{4} y z - p_{4} y - p_{4} z + p_{4} - p_{5} y z + p_{5} z - p_{6} y z + p_{6} y + p_{7} y z$$"
      ],
      "text/plain": [
       "-p₀⋅y⋅z + p₀⋅y + p₀⋅z - p₀ + p₁⋅y⋅z - p₁⋅z + p₂⋅y⋅z - p₂⋅y - p₃⋅y⋅z + p₄⋅y⋅z -\n",
       " p₄⋅y - p₄⋅z + p₄ - p₅⋅y⋅z + p₅⋅z - p₆⋅y⋅z + p₆⋅y + p₇⋅y⋅z"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phix = [1-x, x]\n",
    "phiy = [1-y, y]\n",
    "phiz = [1-z, z]\n",
    "phi = []\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        for k in range(2):\n",
    "            phi.append([phix[i]*phiy[j]*phiz[k]])\n",
    "phi = Matrix(phi)\n",
    "X = Matrix([p0, p1, p2, p3, p4, p5, p6, p7])\n",
    "ref = (Transpose(X)*phi)[0, 0]\n",
    "dx = expand(diff(ref, x))\n",
    "dy = expand(diff(ref, y))\n",
    "dz = expand(diff(ref, z))\n",
    "dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5IAAAARBAMAAABdk47KAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMkSJuyJ2q1Rm\n75mlXrWwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGX0lEQVRYCe1WX4hUVRz+7s7MdXZn171FbCmY\nMiaRgYghFCVuBf2jcqEULGjXepCNWNe3eqj2IemhaMcCYYW2EV/0bcSXLIqBQAu0TAx6HF80Qpad\nULQMt9/5nXPu+XdPrKxvdR7u3PN93/nO9/vdP3OB/8d/pAMnInXG8EC+aKG1stawJks8vVWvSN4I\nvMRw0eWx1DE8amSIWAUx3KxUZ4sWWiuXkNdykae36hXJG4GD7W4TEEsdwxexbayCGB5YLlporVxC\nXstFnt6qVyRvBA62u01ALHUMt7bd8fazBzOeJzNHdo4DL/1SbtOcKwjIANdydjDyqIFy9+UK5ryV\n3R/MHGIBVhzcd6CBypaR5+UcCEhFWLiS+16h0HX3CtaFebDbJPY0VStS7RTgdu0mr+vop5YZgWLc\n7gzGDrfSSd57CN81NiF5Z/aSmHIFPhniWs4ORh4KXXdPrl04b9r3KHazIBnrf6S0BnfPzMmEBAYk\nC21cyz0vpQsNtLtXsI7kwf9ahiT1TqYbBYtMHS7ppVYdUFfSLCroTLW1Hv3zveNjwC7syb5Egiq3\nTVQQkIR5uJaL+JbcEeIJmrruRs6kduE6Pupp4idMvJchzdJubQoZPhR6HgEZ4FrueSkdpMH474G7\nKSydpVeCjmRggBoTLcOQaifZjeeeOV20SMVIx9ueo5VakKoD6kpy9vTYHW72Ox8S474EN5B2T+Ir\nUAfu5xBPZ3i9Xt9Sr68PyBAHSB56OcIXD5Ov5e7ImWSXtF5fe65eb2aDLewZmupbQ93s4buqvBfJ\nWhF307BPhjhA8gIv2+De9vJJ290t+DIuiUb4fUCtGy/DIp3yRheaziInb7auetZujJOaSdkBB1+5\nsHDNyS6yilHpYtnZzRjM6PxPBpriKN+uARnilaZQy2HLLQMcZ9p2Z4Dlkqw0GeE7EqMZ3i1Nla8T\nNNgQ+CtMykNAKs7CpTzwsoTLx/rnaeq4i4dPduMivRLotEkHu14M0ZV0myTmco1NMqrwVx/nWVA7\nF9nbTER9NmlSK1JmlP+TvKhv+/a9tMjJznukG7H6whUMtOheviIQejrF4MICMsSVnNdYckvYkRfL\ncc/lilQuso43kNwsjVSvkmjbMB1qTTro4ZMhruSBl1IKA4CeSc89L7gDfib9PnSwX1zJWBkWqTYS\n3ehs50m4CCKGaLrnaFJLUnZAvV15UQ21EVoVdqangzN3XUepjWO989iPWgd9wp8L88kQ13KxBEZu\nCzO+kq57LpekdpF1bED/z5S9i1r2BXozfAs0eYE4+KQmclzLPS+tUwarhn33vOAMyXna3u9DVj5G\nVzJWhkXqnUQ3stdmLhQu4hiDJ2Yzz9GkZlJ1QF9JWeQQTQs6MzhZPl/uotSo3Fg2X23jyZM/bBZZ\nuDCfDPBczvFzuWMAcSU991zOZO4i67iGF4ZBiTDaeh8fo/zYW5+Ju1AOj9QwclzLPS9X2PsAfPe8\nYODTBvJIBk6TbrwMi9Q7iW5gH04VLuK8053KRq8xJrUgVUZyNDhA32NFnRndOSv+kgYaydH9s6eB\nTv85bhtX4JNk6eK5nOPncjHLhXyxPPdcLq+k3pTzVq7O/Eb8BLByZsXlNpbjzDrWi4NPasLgWu55\neUJ6u3ruVl7xdtWRDPw9XcloGRapdxLdoPNtU0WLOO/qFq56pEktSJWRXAyO8mRhZ+jVS1/cVzBo\n7noCaHAFBWQMF2sceS6UF0vw3hByl+S89P9Co78jjsGIkSEe8ZJCVG/41lZeTLdyNoeTBl3JosFV\nF5Bc3iGs9jvLHhyDmv6352hSu6TBkXa8NWq6Qfyqb1dbcY+YFJAxPJAbofzisc3lObs7ZDJFDP2/\n0NiBl8WPP2JkiEe8WLhquOq3EHle+pecbuf75nBtbu6mgXNeNqmI5PJ+xXRmafNTjtHTpsfOHSa1\nSxocPWPuEjVL+N6c1V+sriZGRvAQlgi+cW2drUOSXiv0Qffg13uLVsXICB7CjJSmeumbKhwy78P4\n3O2+KqNa+EzGSInvEv8SBYNjVCZrawo4xEjGMdAuWoPjfzQITy/SF1Y4YmQED2GJfLL1aOhNr1Xe\nOiBrW+lrBAMLC/MFi2JkBA9hiZTH35wqcFeRnjpwxCVlUpz6S9xj/oiREk8nxvwFYi5j4McJ95aR\n0hip8NKI4/gPRH6tuLQ/tQwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- p_{0} x z + p_{0} x + p_{0} z - p_{0} + p_{1} x z - p_{1} z + p_{2} x z - p_{2} x - p_{2} z + p_{2} - p_{3} x z + p_{3} z + p_{4} x z - p_{4} x - p_{5} x z - p_{6} x z + p_{6} x + p_{7} x z$$"
      ],
      "text/plain": [
       "-p₀⋅x⋅z + p₀⋅x + p₀⋅z - p₀ + p₁⋅x⋅z - p₁⋅z + p₂⋅x⋅z - p₂⋅x - p₂⋅z + p₂ - p₃⋅x⋅\n",
       "z + p₃⋅z + p₄⋅x⋅z - p₄⋅x - p₅⋅x⋅z - p₆⋅x⋅z + p₆⋅x + p₇⋅x⋅z"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5YAAAARBAMAAABUeC6wAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMkSJuyJ2q1Rm\n75mlXrWwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHU0lEQVRYCe1XS2xUZRT+pp1Op0xLr4YQS8Q2\nxcSIxhAMCUQMlYXxEaWJgiEmdpBErIsCKx+JMgtZITK4ICmRUiAxEDdFFj4IOomJuPABhoXEhQUN\nKCHQGgiPGuo553//dy7vnfwJt/f/zjnf+c45d+5/Ae6s/00HvsyoNAtPuZeqKejmgXTWNHLz7Coy\ng/G21nHVTFcv4BZ0ZFSGLDwl5BZyp7jqZL1uHWmyDCSD8bbWoVJnZMrQpeFb0JGVLwtPCbmF3Cmu\nO7MErqOfy956emsizcsN7l7ej46/UTxMe5lZbEzjxl0Ymla9P7hD7ugiuT2kY+v6LdXcUBUbVTpy\ncfRiNJEuyLpHasgzQkLxwlSHhnEHI5QUMZrKbqgOJUPS0yWqL6uarKCwpkhHaFQZy7tGCmvkdjq+\nrc5F+QBaxmgvlcXGNG7chaEw5TGskju6SG6H5MqtCxq7CqVD2G48YOmV0eAuyLpHasgzQkLxwlSH\nhnEHR5IiRlPZjdQBJUPS0yWqL6uajCBNlqGjTsHFkdloHWvpLwMrsDr5ojhyHM08W64sNhIU4dad\nTdjQMIyfsWWLbKQHCsl3UQuTwnipsmFqF94UM10UPRYbo8FdkHV3atidl0MKJ4/F4sVD0Qy8l8DS\nMK7g/lNJLMljHNpB4nQjvDrQbKQGjLoOUNukh5KdLgp/5qkfsjOxsV4QXvs4rsnpKOy9Kwmmdfej\nvO7PYQKF8e/wFZDgASCHo2ivvtrd/Xh39+zYiBSu3HOzmGtuT9I+gtX3lkvDKHR3z/qlu3tYkGmn\n5jBxA9WatFVxCUHu53cZY0AjQco9yMruoY7XsTMQH9BMr0zpUlkD+L7a1DWBpCDHaZzUjQjqSPBN\nEgQFdaA0rmWE9fVNDgdBQSY2qsaHQQ29G4OaAh0zJicv+NMyzw6axtF8aB7aE0Iu0T/aL+mhv/yU\npoxp3LizN9CX4N3GXlBR5qxmBCWaJegRoUtfkj/HG1mSG/v4XowKVTQqyLg7Ncrd1zcfSyoU6YtX\nRKKmkr8YZWV4arl1jJwCSS7HCXq72EbI70EqQ241dykIkkyqjulSNsswS/CXn+BtEOQyKaNoD4K2\nooP3TGY67HRMWbp0bWRkZ3pnzEHnsXNoGwGky/QS3Ma45IuNady4cwiwErkrnT34je9Vbkb0LOUR\nOY6WQ2yVxfSjajhiNLgLMu5WjXb39I3Ow6beSLwiYprG3uJ5wNAILuyg3yXUU2uMXg76XcJU5tXR\n0pnEQcKo6tjMs/SeVNXb0aXiEhTvMiljOkg6GA7E6SihVL/ghlH8OO0iGmvYS988m9E8jHc4u+SL\njWncuItgPILWw/QjWcg7lZsRmWUp2Y6WhLgpwiymT3iW2mhwG2TdrRpxZz8P4VdfKF4RCY28+oKs\nCp7ZY7Iao8eYOwLbCK+OZ2mWmXXk99IslQxTh9T3yuCxOMhlYmO9oAmHmw57OuirqG7B7WvyR/Lj\naKw2TTSPFWv0vDb9a3sVG00PLW7dlfwLeK5nSYIFvFO5GZFZ9o2swwegB76vpnzpyjTgWWqjMdgg\n6y7VW3f28xDKF4lXREJDlUVZBW550GY1kjzGjyjGNMKro0yzzKyjkBs3MkwdIng9DsZBLhMbtXY/\nqHgFLxi8jg7siI0quG/5EB8pbdXcns1D9FGV+2SIDhjdq9iYwq27kDWdH/yLD6+jvJMeCCKznDHY\ncboGzOg/W2GrLKaXWWqjhl2QdZfqrTv7eUhhFJF4IVI0GKCNn1XD9I6NJHmM/I41jXB1tFZollGQ\nZBJh39MstQwB6aIEU0eiID/TkkqdoNwE2nujmpwO5Om7LZyWyriS/uTOUagRwAcFLcmXNqZx5S7R\ndG4AdF7SCaVnKYg+L8WFzy67mF5maRG5CYKUu1Vj3T3kJy88VoPWUbG6rIodxQkXlcqBTSNkVVzS\nQwnqAM3SLsdInwnUwyrNMlpS3w50ut7GmQKjipag82ivWbJYB+jhtcsrmI44Wvo7lm936tT38CY2\nEhThxp29QecDYL5jcxWDeLOc0pWTQZOJltCbD1MF8VVoJMi6u6zG3SFN5aYeHZxWg2V4EZbGss/s\nKcpJQoA1WkY6LTfVbCNcHV+fOXtaJ7JBsuc6SmfOXKlpq/kj9f2KTfoJsEE2E5zRxKimPOx+XKYm\npwMNZeNtjLKnXzOtIf7/pVoL8ba5RcqoLR7uu6OTn+aW0dKwZRAEpUMGaKw1V829pscBC+gbF3QN\nd6Vj8f7PdbPoo8uJV2ryD+1fSx92XlZhb6y0HNbZAqNinI9tSR0u2KmEjCoIxfh3qfAV8pbnZHUy\nOaNWo5vyEj4zQLomtNXqGvf9w70tnKDPKbWWD/WYW6SM2uLhvntpEX1QAAMDprVQSPGPy2aATYN/\nWnZN/+GiPQ7iOy/oGu5Kx++Tk4YgraZtcnIMPo1iz/e/UdFBvlFLenLLbiDNhdZ1n2YH4eDlEaND\n/dVtGihnBxWs0YSqoNb6A9EdplefXiLyP5lu/JXKB0/6AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- p_{0} x y + p_{0} x + p_{0} y - p_{0} + p_{1} x y - p_{1} x - p_{1} y + p_{1} + p_{2} x y - p_{2} y - p_{3} x y + p_{3} y + p_{4} x y - p_{4} x - p_{5} x y + p_{5} x - p_{6} x y + p_{7} x y$$"
      ],
      "text/plain": [
       "-p₀⋅x⋅y + p₀⋅x + p₀⋅y - p₀ + p₁⋅x⋅y - p₁⋅x - p₁⋅y + p₁ + p₂⋅x⋅y - p₂⋅y - p₃⋅x⋅\n",
       "y + p₃⋅y + p₄⋅x⋅y - p₄⋅x - p₅⋅x⋅y + p₅⋅x - p₆⋅x⋅y + p₇⋅x⋅y"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(x_1 - x_0 + x_2 - x_3 + x_4 - x_5 + x_7 - x_6)\\eta\\beta + (x_0 - x_2 + x_6 - x_4)\\eta + (x_0 - x_1 + x_5 - x_4)\\beta + (x_4 - x_0)\\\\\n",
    "(x_1 - x_0 + x_2 - x_3 + x_4 - x_5 + x_7 - x_6)\\beta\\xi + (x_0 - x_1 + x_3 - x_2)\\beta + (x_0 - x_2 + x_6 - x_4)\\xi  + (x_2 - x_0)\\\\\n",
    "(x_1 - x_0 + x_2 - x_3 + x_4 - x_5 + x_7 - x_6)\\xi\\eta + (x_0 - x_1 + x_5 - x_4)\\xi + (x_0 - x_1 + x_3 - x_2)\\eta + (x_1 - x_0)  \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "\\frac{\\mathrm d x}{\\mathrm d \\xi} & \\frac{\\mathrm d x}{\\mathrm d\\eta} & \\frac{\\mathrm d x}{\\mathrm d\\beta}\\\\\n",
    "\\frac{\\mathrm d y}{\\mathrm d \\xi} & \\frac{\\mathrm d y}{\\mathrm d\\eta} & \\frac{\\mathrm d y}{\\mathrm d\\beta}\\\\\n",
    "\\frac{\\mathrm d z}{\\mathrm d \\xi} & \\frac{\\mathrm d z}{\\mathrm d\\eta} & \\frac{\\mathrm d z}{\\mathrm d\\beta}\\\\\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0 1 2]\n",
      "  [3 4 5]\n",
      "  [6 7 8]]\n",
      "\n",
      " [[2 5 8]\n",
      "  [1 4 7]\n",
      "  [0 3 6]]\n",
      "\n",
      " [[8 7 6]\n",
      "  [5 4 3]\n",
      "  [2 1 0]]\n",
      "\n",
      " [[6 3 0]\n",
      "  [7 4 1]\n",
      "  [8 5 2]]\n",
      "\n",
      " [[0 3 6]\n",
      "  [1 4 7]\n",
      "  [2 5 8]]\n",
      "\n",
      " [[6 7 8]\n",
      "  [3 4 5]\n",
      "  [0 1 2]]\n",
      "\n",
      " [[8 5 2]\n",
      "  [7 4 1]\n",
      "  [6 3 0]]\n",
      "\n",
      " [[2 1 0]\n",
      "  [5 4 3]\n",
      "  [8 7 6]]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "idx = np.array([\n",
    "        [0, 1, 2, 3],\n",
    "        [1, 2, 3, 0],\n",
    "        [2, 3, 0, 1],\n",
    "        [3, 0, 1, 2],\n",
    "        [0, 3, 2, 1],\n",
    "        [1, 0, 3, 2],\n",
    "        [2, 1, 0, 3],\n",
    "        [3, 2, 1, 0],   \n",
    "    ])\n",
    "\n",
    "idx0 = np.array([\n",
    "        [0, 1, 2, 3, 4, 5, 6, 7, 8],\n",
    "        [2, 5, 8, 1, 4, 7, 0, 3, 6],\n",
    "        [8, 7, 6, 5, 4, 3, 2, 1, 0],\n",
    "        [6, 3, 0, 7, 4, 1, 8, 5, 2],\n",
    "        [0, 3, 6, 1, 4, 7, 2, 5, 8],\n",
    "        [6, 7, 8, 3, 4, 5, 0, 1, 2],\n",
    "        [8, 5, 2, 7, 4, 1, 6, 3, 0],\n",
    "        [2, 1, 0, 5, 4, 3, 8, 7, 6],\n",
    "    ])\n",
    "print(idx0.reshape(8, 3, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0.  1.  2.]\n",
      "  [ 3.  4.  5.]\n",
      "  [ 6.  7.  8.]]\n",
      "\n",
      " [[ 8.  7.  6.]\n",
      "  [ 5.  4.  3.]\n",
      "  [ 2.  1.  0.]]\n",
      "\n",
      " [[ 6.  3.  0.]\n",
      "  [ 7.  4.  1.]\n",
      "  [ 8.  5.  2.]]\n",
      "\n",
      " [[ 2.  5.  8.]\n",
      "  [ 1.  4.  7.]\n",
      "  [ 0.  3.  6.]]\n",
      "\n",
      " [[ 0.  1.  2.]\n",
      "  [ 3.  4.  5.]\n",
      "  [ 6.  7.  8.]]\n",
      "\n",
      " [[ 8.  7.  6.]\n",
      "  [ 5.  4.  3.]\n",
      "  [ 2.  1.  0.]]\n",
      "\n",
      " [[ 0.  1.  2.]\n",
      "  [ 3.  4.  5.]\n",
      "  [ 6.  7.  8.]]\n",
      "\n",
      " [[ 2.  5.  8.]\n",
      "  [ 1.  4.  7.]\n",
      "  [ 0.  3.  6.]]]\n"
     ]
    }
   ],
   "source": [
    "a = np.array([\n",
    "        [0, 1, 2],\n",
    "        [3, 4, 5],\n",
    "        [6, 7, 8]])\n",
    "\n",
    "I = np.array([0, 2, 3, 1, 0, 2, 0, 1])\n",
    "b = np.zeros((len(I), 3, 3))\n",
    "for i, k in enumerate(I):\n",
    "    b[i] = np.rot90(a, k=k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0. ,  0.5,  1. ])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linspace(0, 1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Hamiltontian "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
